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ABSTRACT 


It  is  shown  that  solutions  to  the  plane,  radial,  steady,  quasi¬ 
static  flow  of  a  Bingham  solid  in  a  converging  channel,  and  to  the  axially 
symmetric,  radial,  steady,  quasi -stati c ,  converging  flow  in  a  circular 
cone  do  not  exist,  at  least  if  the  material  is  incompressible,  and  there 
are  no  body  forces . 

The  velocity  and  stress  fields  for  the  corresponding  power  law 
pseudo-plastic  flow  problems  are  obtained.  A  new  numerical  technique  is 
developed  to  solve  the  non-linear  boundary  value  problem  that  arises  in 
the  determination  of  these  fields. 

The  results  are  applied  to  hydraulic  extrusions  and  wire  drawing 


processes . 


. 

jv  |  nf:r)qi'.moD  eril  rio'  et  r^1  ■  -  v 


ACKNOWLEDGMENTS 


The  author  wishes  to  extend  his  appreciation  to  Dr.  J.B.  Haddow 
for  his  guidance  and  supervision  of  this  thesis.  He  also  gratefully 
acknowledges  the  financial  support  for  this  research  by  the  National 
Research  Council  and  by  the  University  of  Alberta. 

Thanks  are  extended  to  the  members  of  the  Mechanical  Engineering 
Shop  who  made  the  experimental  apparatus  and  particularly  to  R.  Haswell 
who  assisted  in  all  stages  of  the  experimental  work. 

Additional  thanks  are  also  extended  to  Dr.  J.M.  Chudobiak  who 
on  numerous  occasions  curtailed  his  own  research  in  order  to  aid  the 
author,  and  to  Miss  H.  Wozniuk  for  her  patience  in  typing  the  thesis. 

The  author  would  like  also  to  express  his  appreciation  to 
his  wife,  Ineke,  for  her  encouragement  and  considerations. 


iv 


' 


TABLE  OF  CONTENTS 


Page 

CHAPTER  I  INTRODUCTION  1 

1 . 1  The  Bingham  Sol  id  4 

1.2  The  Pseudo-Plastic  Material  8 

1.3  Aim  of  the  Thesis  10 

1.4  Consideration  of  Related  Problems  11 

Plane  Viscous  Flow  -  Hamel's  Problem  12 

Axially  Symmetric  Analog  to  Hamel's  Problem  12 

Plane  Rigid  Perfectly  Plastic  Flow  13 

Axially  Symmetric  Rigid  Perfectly  Plastic 

Flow  13 

CHAPTER  II  VISCO-PLASTIC  FLOW  15 

2.1  Plane  Visco-Plastic  Flow  15 

2.2  Axially  Symmetric  Visco-Plastic  Flow  24 

CHAPTER  III  PSEUDO-PLASTIC  FLOW  33 

3.1  Plane  Pseudo-Plastic  Flow  33 

3.2  Axially  Symmetric  Pseudo-Plastic  Flow  42 

3.3  Integration  Procedure  50 

3.4  On  the  Inclusion  of  Inertia  Effects  52 

3.5  On  the  Neglect  of  Inertia  Effects  54 


v 


, 


"/  •I  r  ,•  -o:w  \b  • 


W.  i  6fc‘  0  b 

' 


TABLE  OF  CONTENTS  (continued) 


Page 

CHAPTER  IV  APPLICATION  TO  WIRE  DRAWING  AND  HYDRAULIC 

EXTRUSION  63 

4.1  Wire  Drawing  68 

4.2  Hydraulic  Extrusion  78 

CHAPTER  V  EXPERIMENTAL  INVESTIGATION  82 

5.1  Experimental  Apparatus  83 

5.2  Results  of  Experiments  88 

The  Specimens  88 

Lead  Extrusions  89 

Wax  Extrusions  90 

5.3  Conclusions  91 

CHAPTER  VI  CONCLUDING  REMARKS  93 

6.1  Results  93 

6.2  Error  Analysis  of  the  Numerical  Results  128 

6.3  Objections  to  the  Power  Law  Model  130 

6.4  Wire  Drawing  and  Hydraulic  Extrusion  132 

6.5  Suggestions  for  Future  Work  141 

BIBLIOGRAPHY  147 


vi 


■■  '  ■  ■  . 


1 '  n  r?v'^  ■ ,  : 


3HIGUJ0HQ0 

nor2UT$x3  d  hr  6-0  <r<  .  £ 


TABLE  OF  CONTENTS  (continued) 


Page 

APPENDIX  NUMERICAL  INTEGRATION  OF  ORDINARY  DIFFERENTIAL 


EQUATIONS 

A.l  Introduction  152 

A. 2  The  First  Order  Ordinary  Differential 

Equation  152 

A. 3  Equations  Governing  the  Constants  and 

Weighting  Factors  154 

A. 4  Determination  of  the  Constants  and 

Weighting  Factors  158 

A. 5  The  System  of  First  Order  Ordinary 

Differential  Equations  165 

A. 6  Gill's  Variation  of  the  Runge-Kutta 

Fourth  Order  Method  169 

A. 7  Brief  Discussion  on  the  Accuracy  of  the 

Method  174 

A. 8  Subroutine  DRKGS  177 


VI  1 


. 


LIST  OF  TABLES 


Table 
A.  1 


Typical  Terms  and  their  Coefficients 


Page 

169 


vi  i  i 


LIST  OF  FIGURES 


Figure  Page 

1.1  Stress-Strain  Rate  Relationships  for  Simple  Shear  4 

2.1  Visco-Plastic  Flow  in  a  Converging  Channel  16 

2.2  Visco-Plastic  Flow  in  a  Circular  Cone  25 

3.1  Pseudo-Plastic  Flow  in  a  Converging  Channel  34 

3.2  Compression  of  a  Pseudo-Plastic  Material  Between 

Inclined  Plates  42 

3.3  Pseudo-Plastic  Flow  in  a  Circular  Cone  43 

3.4  Pseudo-Plastic  Flow  in  a  Circular  Tube  56 

3.5  Velocity  Components  in  a  Circular  Cone  of  Small 

Semi -Angle  a  60 

4.1  Flow  Curves  for  the  Power  Law  Model  64 

4.2  Flow  Curve  for  the  "Rigid"  Pseudo-Plastic  Material  66 

4.3  Cross-Section  of  a  Conical  Reducing  Die  67 

5.1  Typical  Half  of  the  Die  85 

5.2  The  Separation  Plate  86 

5.3  The  Piston  and  the  Piston  Housing  87 

6.1  Possible  Negative  Normal  Stresses  at  the  Walls  134 

6.2  The  Cone  of  Finite  Length  142 


ix 


9boM  WfiJ  "OWOl  9  rtf  icl  29V'*jO  won 

9ra  gnr '  sP  Tso  noD  s  -o  d 

sH  c 

gnrauoH  notsr9  aria  bno  no.  3  9:1 
a T  T bW  9fi^  ^6  zezzo'ilZ  f  rm o  1  :k-  -j  5r,  ,' 


LIST  OF  FIGURES  (continued) 


Figure  Page 

A.l  Geometrical  Construction  of  6y  154 

Results  of  Plane  Pseudo-Plastic  Flow  in  a 
Converging  Channel 

6.11- 1  m1 ' (0)  versus  n  99 

6.11- 2  m( 0)  versus  0:  of  15°  100 

6.11- 3  m(0)  versus  0:  of  30°  101 

6.11- 4  m(0)  versus  0:  of  45°  102 

6.11- 5  m ( 0 )  versus  0:  of  60°  103 

6.11- 6  t  Q(a)  r2n/yQn  versus  n  104 

6.11- 7  Tr0(a)  r2n/uQn  versus  a:  a  =  0°-30°  105 

6 . 1 1 -  8  r  Q(a)  r^n/yQn  versus  a:  a  =  30°-60°  106 

ru 

6.11- 9  aTr0(a)  r2n/pQn  versus  a:  a  =  0°-25°  107 

6.11- 10  a..r^n/yQn  versus  0:  a  =  15°,  n  =  0.1  108 

*  vJ 

6. 11- 11  a. .r2n/uQn  versus  6:  a  =  30°,  n  =  0.1  109 

*  J 

6. 11- 12  a.  .r^n/yQn  versus  0:  a  =  45°,  n  =  0.1  110 

*  J 

6. 11- 13  a. .r^n/yQn  versus  0:  a  =  60°,  n  =  0.1  111 

'  J 

Results  of  Axially  Symmetric  Pseudo-Plastic 
Flow  in  a  Circular  Cone 

6.111- 1  m ' ' ( 0 )  versus  n  113 

6.111- 2  m( 4>)  versus  <J> :  a  =  15°  114 

6.111- 3  m( 4>)  versus  <j):  a  =  30°  115 

6.111- 4  m( 4>)  versus  <J>:  a  =  45°  116 


x 


■ 


r.A 


LIST  OF  FIGURES  (continued) 


Figure 

6.111- 5  m(cf>)  versus  4>:  a  =  60° 

6.111- 6  Trcj)(a)  r3n/yQn  versus  n 

6.111- 7  T^(a)  r3n/yQn  versus  a:  a  =  0°-30° 

6.111- 8  T^(a)  r3n/yQn  versus  a:  a  =  30°-60° 

6.111- 9  ai^(a)  r3n/yQn  versus  a:  a  =  0°-20° 

6.111- 10  a..  r3n/yQn  versus  :  a  =  15°,  n  =  0.1 

*  J 

6.111- 11  a..  r3n/yQn  versus  cj>:  a  =  30°  ,  n  =  0.1 

*  J 

6 . 1 1 1  - 1 2  a..  r3n/yQn  versus  4> :  a  =  45°,  n  =  0.1 

*  J 

6.111- 13  a.  j  r3n/yQn  versus  <(>:  a  =  60°  ,  n  =  0.1 

*  J 

6.111- 14  m((|))  versus  4>:  a  =  5°,  n  =  0.1 

6 . 1 1 1  - 1 5  a..  r3n/yQn  versus  (£:  a  =  5°,  n  =  0.1 

*  J 

Results  of  Wire  Drawing  and  Hydraulic  Extrusion 
6.A-1  Maximum  Percentage  Reduction  versus  a 
6.A-2  Drawing  Stress  x  D-|3n/yQn  versus  a 
6 -A- 3  Pram  D13n/uQn  versus  a 


Page 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 
127 

138 

139 

140 


xi 


NOMENCLATURE 


Qij 
Sij 
-  P 

d.  . 
1J 


f . 


Vi 


y 

k 

Y 

I 

J 

D 

Dt 


n 

a 

D- 


components  of  a  stress  tensor 
components  of  a  deviatoric  stress  tensor 
hydrostatic  part  of  the  stress  tensor 
components  of  a  strain  rate  tensor 
densi ty 

body  forces  per  unit  mass 
components  of  a  velocity  vector 
components  of  an  acceleration  vector 
(apparent)  coefficient  of  viscosity 
yield  stress  in  pure  shear 
tensile  yield  stress 
/icT  .d. .  a  strain  rate  invariant 

*  U  *  U 

/i  s  . . s  .  a  stress  invariant 

■  vj  *  J 

material  derivative 
covariant  differentiation 
parameter  in  power  law  model 
semi-angle  of  the  channel  or  cone 

exit  diameter  of  the  circular  cone  in  wire  drawing  and 

hydraulic  extrusion 

volume  flow  (per  unit  length) 


2n ,  ^n 


a..r  /yQ  components  of  a  non-dimensional  stress  tensor  for  plane  flow 
•  J 


a..r'jn/yQn  components  of  a  non-dimensional  stress  tensor  for  axially 
^  J 

symmetric  flow 


XT  1 


*  ozrr  \  z  v  t  ft  <  i.  »  ,.nc  o 

'<U2f  ‘-VI  ,flt2  '  .03 

2?6i ,  o  •  n.  k  jc  osmcrt  ybod 


/  no  JMfsPj  fc  ^  ;.:  ) 


nor.:  . r  >  9r<  tl  b  Jnsnj  ;oo 


NOMENCLATURE  (continued) 

P  Dn  /uQn  non-dimensional  ram  pressure  in  hydraulic  extrusion 
ram  1  r  r  J 

Drawing  stress  x  D-j  /yQn  non-dimensional  drawing  stress  in  wire 
drawing 


XT  1  1 


CHAPTER  I 


INTRODUCTION 

Analysis  of  the  deformation  and  flow  of  matter  is  facilitated 
by  the  assumption  that  matter  is  a  continuum  and  by  further  introducing 
the  broad  classification  of  matter  as  being  either  a  fluid  or  a  solid. 

A  continuum  is  a  hypothetical  representation  of  matter  as  a 
continuous  medium;  a  study  based  upon  the  assumption  of  a  continuum  is 
termed  a  phenomenological  approach.  A  fluid,  according  to  the  classical 
definition,  is  a  substance  that  flows  under  the  action  of  any  anisotropic 
stress  system,  no  matter  how  small  (1)*.  A  material  is  said  to  flow  if  it 
deforms  continuously  with  time  (1).  A  solid,  according  to  the  classical 
definition,  is  a  substance  that  requires  a  definite  non-zero  anisotropic 
stress  level  to  produce  flow;  a  stress  below  this  level  produces  defor¬ 
mation,  which  may  be  reversible  or  irreversible,  but  no  flow. 

The  classification  of  a  "real"  substance  as  a  solid  is  to  some 
extent  an  idealization;  this  becomes  more  evident  as  experiments,  employing 
more  refined  measuring  techniques,  show  that  the  application  of  an  aniso¬ 
tropic  stress  system  to  a  "real"  solid  results  in  continuous  deformation. 
Indeed  the  motto  of  the  Society  of  Rheology  is  "Everything  flows"**  and 
this  is  termed  by  Reiner  (1)  as  "the  second  axiom  of  rheology".  The  state- 

*Numbers  without  decimal  enclosed  by  brackets  designate  the  references 
listed  in  the  Bibliography. 

'f 

**The  motto  of  the  Society  of  Rheology  had  its  origin  in  the  strictly 
philosophical  deduction  by  Hericlitus  (3),  that  "Everything  flows". 
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ment  "Everything  flows"  is  justified  since  the  essential  difference  between 
the  tendency  to  flow  of  a  fluid  and  a  solid  is  quantitative  and  not  quali¬ 
tative.  Thus,  according  to  Reiner  (2),  it  is  equally  as  justifiable  to 
claim  that  "Everything  is  solid". 


The  analysis  of  the  deformation  and  flow,  resulting  from  the 


application  of  forces  to  a  material  in  bulk,  is  based  on  the  following 
equations : 


(1.1) 


(1.2) 


(1.3) 


and  equations  resulting  from  the  first  and  second  law  of  thermodynamics, 
which  are  not  included  since  thermodynamic  effects  are  neglected  in  the 
analysis  presented  in  this  thesis.  Equation  (1.2)  is  valid  since  the  body 
moments  per  unit  mass  and  couple  stresses  are  assumed  to  be  zero.  The 
relations  (1.1),  (1.2),  and  (1.3)  provide  seven  equations  for  the  thirteen 
unknowns  a..,  v. ,  and  p,  if  the  distribution  of  the  body  force  per  unit 

I  J  I 

mass,  f. ,  is  known.  The  problem,  therefore,  remains  six  degrees  inde¬ 
terminate.  To  proceed  further,  the  material  response  has  to  be  considered; 
consequently  what  are  known  as  constitutive  equations  are  introduced. 
Eringen  (4)  notes  that  "The  character  of  the  material  is  brought  into  the 
formulation  through  appropriate  constitutive  equations  for  each  material 


. 
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with  the  constitutive  variables  being  restricted  in  their  regions  of  defi- 
ni tions . " 

Constitutive  equations  cannot  be  formulated  arbitrarily,  but 
must  satisfy  certain  physical  and  mathematical  requirements  based  on  the 
axioms  of  constitutive  theory  as  described  by  Eringen  (4).  The  constitutive 
equations  are  relations  between  the  stress  and  strain  tensors  and  their 
time  derivatives.  For  an  arbitrary  element  these  equations  must  be  inde¬ 
pendent  of  any  superimposed  rigid  body  motion.  In  the  constitutive  equations 
certain  rheological  properties  appear  as  parameters;  these  are  elasticity, 
viscosity,  and  plasticity.  These  properties  are  termed  by  Reiner  (5)  as 
fundamental  and  other  properties,  which  are  combinations  of  the  fundamental 
properties,  as  complex  properties.  Some  examples  of  complex  properties 
are  retardation,  relaxation,  delayed  elasticity,  visco-plasticity,  visco¬ 
elasticity,  and  anelastici ty . 

A  "real"  material  exhibits  all  of  the  fundamental  rheological 
properties  to  some  extent  (6).  The  inclusion  of  certain  specific  rheological 
properties  in  the  constitutive  equation  of  a  "real"  material  is  justified 
if  it  is  observed  experimentally  that  these  properties  dominate.  This 
thesis  considers  visco-plastic  and  pseudo-plastic  materials.  These  are 
inelastic  materials  for  which  plastic  and  viscous  properties  dominate. 
According  to  the  classical  definitions,  the  visco-plastic  material  con¬ 
sidered  in  this  thesis  is  an  example  of  a  solid,  and  the  pseudo-plastic 
material  is  an  example  of  a  fluid.  The  relationships  between  shearing  stress, 

t  ,  and  shearinq  strain  rate,  d  ,  in  simple  shear  for  these  materials  are 
xy  3  xy 

indicated  in  FIG.  1.1,  where,  for  the  purpose  of  comparison,  the  Newtonian 
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viscous  liquid  is  also  shown. 

Shear  Stress 


FIGURE  1.1 


Stress-Strain  Rate  Relationships  for  Simple  Shear 


The  remainder  of  this  chapter  presents  a  brief  discussion  of  the  Bingham 
solid  and  the  pseudo-plastic  material,  respectively. 

1 . 1  The  Bingham  Solid 

In  FIG.  1.1  the  flow  curve  for  a  Newtonian  viscous  liquid  is 
shown.  The  Newtonian  liquid  is  a  viscous  material  and  its  constitutive 
equation  in  simple  shear  is 


T 

xy 


n  "'fit  ij  w  y  t:  wol-f  ant  r  r)  r  nl 
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The  coefficient  y,  which  depends  only  on  temperature  and  pressure, and  not 
on  rate  of  deformation,  is  called  the  Newtonian  coefficient  of  viscosity, 
and  it  completely  characterizes  the  fluid  if  the  fluid  can  be  considered 
incompressible. 

The  fundamental  rheological  properties  may  be  associated  with 
the  manner  in  which  mechanical  energy  is  dissipated  by  the  material.  A 
material  which  is  classified  as  perfectly  elastic  does  not  dissipate  energy*; 
the  Newtonian  viscous  liquid  dissipates  energy  by  the  transfer  of  momentum 
on  a  molecular  scale.  Gases  and  liquids  of  low  molecular  weight  are  con¬ 
tained  in  the  viscous  category  (7). 

An  inelastic  dissipative  material,  whose  behavior  is  not  even 
approximately  described  by  the  Newtonian  constitutive  equation,  is  paint  (8). 
Paint  is  required  to  possess  the  following  practical  properties: 

(1)  it  must  brush  without  too  much  effort, 

(ii)  it  must  flow  in  order  that  brushmarks  disappear,  i.e.  it  must 
possess  "leveling"  properties,  and 
(iii)  it  must  not  run  when  applied  to  a  vertical  wall. 

Reiner  (8)  points  out  that  a  low  viscosity  is  required  for  the  first  two 
properties,  and  a  high  viscosity  for  the  third;  the  viscosity  requirements 
seem  to  be  in  conflict.  The  solution  to  this  problem  was  found  by  Bingham 
and  Green  (9)  and  was  presented  by  them  in  a  paper  entitled:  "Paint,  a 
Plastic  Material  and  not  a  Viscous  Liquid".  The  basis  of  the  solution  rested 
upon  the  experimentally  determined  fact  that  the  viscosity  of  paint  is  not 
a  constant  for  different  rates  of  flow. 


*At  least  if  thermal  effects  are  neglected. 
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The  one-dimensional  model,  which  Bingham  and  Green  (9)  suggested 
to  describe  the  behavior  of  paint,  has  as  its  rheological  equation  in 
simple  shear  (10) 


2ydxy  +  T0  sgn 


if 


>  T 


0* 


d 


xy 


=  o 


if 


xy 


<  T 


This  two-parameter  model  is  now  termed  the  Bingham  solid.  A  generalization 
of  this  model  to  three  dimensions  was  suggested  by  Hohenemser  and  Prager  (11). 
In  three  dimensions  Prager  (12)  has  characterized  the  mechanical  behavior 
as  follows:  "Let  a  Newtonian  viscous  liquid,  a  perfectly  plastic  Mises 
solid,  and  a  visco-plastic  Bingham  solid  be  subjected  to  the  same  velocity 
strain.  The  stress  in  the  Bingham  solid  is  then  obtained  by  adding  the 
stresses  in  the  Newtonian  liquid  and  the  Mises  solid".  The  constitutive 
equations  for  this  generalization  are 


s.j  =  2(u.+  y)  dij  .  if  I  *  0. 


(1.4) 


where 


I  =  (2d  d  ) 
mn  mn7 


1/2 


Clearly  the  Newtonian  viscous  liquid  and  the  perfectly  plastic  von  Mises 
solid  are  special  cases  of  the  Bingham  solid;  this  is  consistent  with  the 
third  axiom  of  rheology,  which  states  that  the  constitutive  equation  of  a 
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more  simple  body  can  be  derived  from  that  of  a  less  simple  body  (13). 

The  constitutive  equation  (1.4)  may  be  inverted  as  follows. 
Squaring  both  sides  of  (1.4)  gives 


where 


if  J  <  k, 
J  >  k, 


J  -  (i  smn  smJ1/2 
2  mn  mn 


Substituting  in  the  original  equation  (1.4)  gives 


2y 


d.  . 
ij 


=  { 


0 

o-j) 


ij 


,  if  J  <  k, 

,  if  J  >  k  . 


Few  non-trivia!  problems  have  been  solved  for  the  Bingham  solid. 
Oldroyd  (14)  considered  some  two-dimensional  boundary  layer  and  rectilinear 
flow  problems.  Carlson  (15)  has  treated  the  axial  compression  of  a  disc. 
Carlson  used  a  piecewise  linearized  yield  criterion  introduced  by  Prager  (16) 
in  1961.  Haddow  (17)  extended  Carlson's  solution  to  include  inertia  effects. 
The  linearization  theory  was  extended  by  Haddow  (18)  and  was  applied  by 
Haddow  and  Hrudey  (19)  to  solve  the  flow  problem  of  a  thin  circular  plate 
subjected  to  a  uniformly  distributed  transverse  load.  Prager  (20)  derived 
extremum  principles  for  the  Bingham  solid;  these  extremum  principles  were 
extended  by  Haddow  and  Luming  (21). 
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1 .2  The  Pseudo-Plastic  Material 

Comparison  of  the  Bingham  constitutive  equation, 


s, =  2(u  +  ■  _ ^  )  d..  , 

/  2d  _  d  1J 


mn  mn 


with  that  of  a  Newtonian  viscous  liquid, 


suggests  that  the  Bingham  solid  can  be  considered  to  be  a  material  with  a 
variable  coefficient  of  viscosity  that  depends  on  the  strain  rate  invari¬ 
ant  I , 


u 


s 


which  can  be  termed  the  "apparent"  viscosity  Other  examples  of  materials 
with  a  variable  coefficient  of  viscosity  are  the  pseudo-plastic  model, 
which  describes  a  material  for  which  viscosity  decreases  with  an  increasing 
shear  rate,  and  the  dilatant  fluid,  for  which  the  viscosity  increases  with 
an  increasing  rate  of  shear  (7).  These  effects  were  originally  noted  in 
colloidal  solutions  and  their  explanation  is  due  to  Ostwald  (22).  Ostwald 
hypothesized  that  a  "structure"  in  the  fluid  altered  with  flow  and  introduced 
the  term  "Strukturvi skosi tat"  or  "Structural  viscosity"  to  denote  this 
phenomenom.  Fluids  exhibiting  this  kind  of  variable  viscosity  have  been 
designated  by  Reiner  (23)  as  "Non-Newtonian  fluids";  this  "Non-Newtonian" 


•  i 
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classification  includes  the  Bingham  solid,  pseudo-plastics  and  dilatant 
fluids*  (7). 


Ostwald  (10)  suggested  a  two-parameter  power  law  model  as  a 


means  of  characterizing  the  pseudo-plastic  or  the  dilatant  fluid;  this 
power  law  is  commonly  written  in  simple  shear  as 


(1.5) 


For  this  constitutive  equation  the  apparent  viscosity  is 

n-1 


Thus,  n  <  1  characterizes  a  pseudo-plastic,  and  n  >  1  characterizes  the 
dilatant  fluid.  Reiner  (24)  has  raised  three  major  objections  to  the  use 
of  the  power  law  and  concludes  that  equation  (1.5)  is  not  a  physical  law 
but  rather  only  an  empirical  formula.  An  obvious  defect  in  this  formula 
is  that  the  physical  dimension  of  p  depends  on  the  value  n. 

Other  empirical  equations,  which  describe  pseudo-plastic  behavior 
in  simple  shear,  are 


Prandtl 


Eyring 


*0ften  the  terminology  "Non-Newtonian"  fluids  includes  a  broader  class  of 
materials  such  as  elastico-viscous  materials  and  models  which  show  a 
cross  viscosity. 
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Powell -Eyring  Txy  =  2Adxy  +  B  sinh'1  (2  C  d  )  , 

Williamson  t  =  2Ad  / ( B  +  2d  )  +  y  d  , 

xy  xy'  v  xy'  xy  * 

where  A,  B,  and  C  are  constants  characterizing  the  particular  fluid  be¬ 
havior  (7) . 

1 .3  Aim  of  the  Thesis 

The  purpose  of  this  thesis  is  to  consider  the  plane  flow  of  an 
incompressible  visco-plastic  and  a  pseudo-plastic  material  in  a  converging 
channel,  and  also  axially  symmetric  converging  flow  in  a  circular  cone. 

The  visco-plastic  material  considered  is  the  Bingham  solid  and  the  pseudo¬ 
plastic  material  is  based  on  the  power  law  model.  The  constitutive  equation 
for  the  Bingham  solid,  which  is  applied,  is  the  Hohenemser-Prager  gener¬ 
alization 


2(y  + 


/~2d 


d 

mn  mn 


— )  d .  .  . 
1  iJ 


The  constitutive  equation  which  is  used  for  the  pseudo-plastic  material, 
is  the  power  law  representation  (25) 


s 


•  • 


U 


n-1 

2 


(1.6) 


The  special  case  of  equation  (1.6) wi th  n  =  0  is  the  perfectly  plastic  von 
Mises  solid,  whose  yield  criterion  is 


-sd  bTutt  ifiFuott'YBq  sitt  a  ns  a.ioo  d^n  mb  «;  *• 

I  1 1 1  'JbBO  SI  W  fl|  1 1  k*)  iohnM 


11 


2  smn 


» 


where  k  is  the  yield  stress  in  pure  shear.  When  n  =  1,  equation  (1.6) 
is  that  for  a  Newtonian  viscous  liquid. 

The  visco-plastic  material  as  typified  by  the  Bingham  solid 
possesses  a  yield  limit;  this  is  to  be  contrasted  with  the  pseudo-plastic 
material  which  does  not  possess  a  yield  limit  for  values  of  n  in  the  range 

0  <  n  £  1 . 

However,  the  limiting  case  n  =  0,  the  perfectly  plastic  von  Mises  solid, 
does  possess  a  yield  limit. 

The  rheological  problems  considered  in  this  thesis  may  be  appli¬ 
cable  to  hydraulic  extrusion  and  wire  drawing  technologies  which  involve 
rate  dependent  materials. 

Calculations  are  performed  for  values  of  n  between  0  and  1; 
thus  it  is  possible  to  study,  both  quantitatively  and  qualitatively,  the 
transition  from  a  Newtonian  viscous  liquid  to  a  perfectly  plastic  von 
Mises  solid. 

1 .4  Consideration  of  Related  Problems 

Various  solutions  have  been  presented  for  the  plane  flow  of 
Newtonian  viscous  liquids  and  rigid  perfectly  plastic  materials  in  a  con¬ 
verging  channel,  and  also  axially  symmetric  converging  flow  in  a  circular 
cone.  These  flow  problems  shall  now  be  considered  briefly. 


. 


bMoz  rrTbrtgnfS  aril  yd  bsmqyd  as 
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Plane  Viscous  Flow  -  Hamel's  Problem 

In  a  classic  paper,  and  using  the  prior  assumption  of  radial 
flow,  Rosenhead  (26)  obtained  closed  form  solutions  for  the  plane,  steady 
flow  of  an  incompressible  Newtonian  viscous  liquid  between  two  inclined 
walls,  separated  by  an  angle  of  2a.  Rosenhead  included  inertia  terms 
and  made  a  thorough  investigation  of  the  change  in  flow  characteristics 
due  to  increasing  Reynold's  number. 

The  mathematical  solution  is  extremely  complex  and  is  obtained 
in  terms  of  Weierstrassian  elliptic  functions.  A  principal  result  of 
Rosenhead 's  investigation  is  that  for  every  pair  of  values  a  and  Reynold's 
number,  the  number  of  mathematically  possible  velocity  profiles  is  in¬ 
finite.  The  profiles  may  or  may  not  be  symmetrical  with  respect  to  the 
central  line  of  the  channel.  Further,  Rosenhead  found  that  both  inflow 
and  outflow  could  occur  simultaneously.  A  solution  not  assuming  radial 
flow  has  not  been  obtained  as  yet. 

Axially  Symmetric  Analog  to  Hamel's  Problem 

It  has  been  shown  by  Hamel  (27)  and  others  that  for  axially 
symmetric  converging  flow  of  an  incompressible  Newtonian  viscous  liquid 
in  a  circular  cone,  without  body  forces,  no  purely  radial  flow  can  exist, 
if  inertia  effects  are  considered. 

The  assumption  of  non-radial  flow  leads  to  an  intractable  problem. 

It  is  interesting  to  note  that  a  radial  flow  solution  is  obtained  upon  the 
neglect  of  inertia  terms.  This  solution  has  been  employed  by  Ramacharyulu  (28) 
as  a  first  order  approximation  to  the  elastico-viscous  flow  problem  in  a 
circular  cone.  The  second  order  solution  yielded  a  non-radial  flow  pattern. 


i  -  f 


b  Mtttoti  ow*  tmrtsd  MtfpK  auosafcv  nfinoJwsK  efdraa^qrncort-  ne  o  wor  t 


e  nf  rasfdoiq  wqtt  auooaW-oardasrs  sHJ  oi  noHsmfxo'iqq6  isbtq  Jaitt  s  af 


13 


Plane  Rigid  Perfectly  Plastic  Flow 

In  1924  Nadai  (29)  considered  the  plane,  steady,  quasi-static 
flow  of  an  incompressible  rigid  perfectly  plastic  solid  in  a  converging 
channel,  of  semi-angle  a,  with  perfectly  rough*  walls.  Nadai  assumed 
that  the  stresses  are  a  function  of  9,  but  not  of  r,  and  this  implies 
that  the  streamlines  of  the  associated  velocity  field  must  be  radial 
straight  lines  directed  through  the  virtual  apex  of  the  channel.  Further¬ 
more,  assuming  that  the  channel  is  very  long,  and  thus  neglecting  end 
effects,  enabled  Nadai  to  obtain  the  stress  field  for  this  flow  problem. 
Hill  (30)  obtained  the  associated  velocity  field  of  the  form 


It  may  be  shown  that  g'(0)  at  0  =  a  is  infinite.  This  derivative  cannot 
be  infinite  in  the  Bingham  solid  since  this  would  imply  that  the  shearing 
stresses  at  the  wall  are  also  infinite.  Hill  (30)  has  extended  Nadai's 
solution  to  cover  the  cases  of  rough  walls  with  a  constant  shearing  traction 
mk  at  the  walls ,  where 


0  <_  m  £  1 . 

Axially  Symmetric  Rigid  Perfectly  Plastic  Flow 

In  1955  Shield  (31)  obtained  the  stress  field  for  the  axially 


*The  coefficient  of  friction  between  a  perfectly  rough  surface  and  the 
deforming  rigid  perfectly  plastic  material  is  large  enough  that  the 
shearing  yield  stress  can  be  developed  at  the  surface. 


rf  ^ngjgnoa  6  fttjw  af  ffei  dgoi m  mbo  9f  t  yevoo  o:  norJufoa 


won  oM<:*n  lloa'.-oq  br?n  onJam-TiZ  y,r fsrxA 


14 


symmetric  converging  flow  of  an  incompressible  rigid  perfectly  plastic 
material  in  a  circular  cone  of  semi-angle  a.  The  wall  friction  is  as¬ 
sumed  to  be  constant  at  the  value  of  mk,  where 

0  <  in  <  1, 

and  k  is  the  yield  stress  in  pure  shear.  The  velocity  field  is  radial, 
and  consequently  of  the  form 


Shield  did  not  obtain  a  closed  form  solution  to  the  governing 
differential  equation,  which  is  in  terms  of  the  shearing  stress  t^,  but 
solved  it  numerically.  Shield  applied  the  results  of  the  stress  field  to 
wire  drawing  in  a  manner  to  be  discussed  in  Chapter  IV. 

As  mentioned  previously  in  this  chapter,  it  is  desired  to  study 
wire  drawing  and  hydraulic  extrusion.  Most  materials  used  in  this  process 
exhibit  both  plastic  and  viscous  effects. 

The  remaining  chapters  are  concerned  with  the  application  of 
the  Bingham  solid  and  the  pseudo-plastic  material  to  the  above  mentioned 
plane  and  axially  symmetric  flow  problems. 
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CHAPTER  II 


VISCO-PLASTIC  FLOW 

This  chapter  considers  plane  flow  of  an  incompressible  homo¬ 
geneous  isotropic  rigid  visco-plastic  material,  the  Bingham  solid,  in  a 
converging  channel,  and  axially  symmetric  converging  flow  in  a  circular 
cone. 

It  will  be  shown  that  the  only  possible  radial  flow  solution 
is  the  trivial  solution,  with  zero  radial  velocity,  unless  it  is  assumed 
that  slip  occurs  at  the  walls  and  axially  symmetric  flow  results  for  the 
plane  problem  and  spherically  symmetric  flow  for  the  axially  symmetric 
problem. 

The  assumption  of  a  more  general  velocity  field,  which  includes 
both  radial  and  tangential  components,  leads  to  extensive  and  cumbersome 
equations,  which  cannot  be  integrated  analytically. 

2.1  Plane  Visco-Plastic  Flow 

The  plane  flow  of  an  incompressible  Bingham  solid  In  a  converging 
channel  is  considered.  Polar  coordinates,  (r,6),  which  are  taken  relative 
to  the  axis  of  symmetry,  are  used  to  describe  the  flow,  and  the  channel 
is  bounded  by  the  lines  6  =  ±  a.  It  is  assumed  that  the  flow  is  steady, 
quasi -static,  and  directed  radially  through  the  virtual  apex  of  the  channel, 
0,  see  FIG.  2.1,  and  that  the  channel  is  long  enough  so  that  the  velocity 
and  stress  fields  are  independent  of  conditions  at  the  ends. 
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Axis  of 
Symmetry 


Visco-Plastic  Flow  in  a  Converging  Channel 


The  components  of  the  velocity  strain  tensor,  defined  by 


are 


d.  . 
ij 


=  2(vi,j 


v .  . ) 

J.V 


(2.1) 


» 


(2.2) 
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re 


2 r  ae 


The  condition  of  incompressibility  is 


which  using  equations  (2.2)  is 


3v 


r 


(2.3) 


Equation  (2.3)  requires  that  the  velocity  field  be  of  the  form 


Substitution  of  equation  (2.4)  in  equations  (2.2)  yields 


(2.5) 


where  1  denotes  differentiation  with  respect  to  0. 

The  components  of  the  stress  tensor,  which  do  not  vanish  identi¬ 
cally,  are  a  ,  ae>  and  xrQ.  Assuming  steady,  quasi-static  flow,  and  neg- 
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lecting  body  forces,  the  equilibrium  equations  are 


do  3t  - 

Vs  K'ft 

r  —  +  —  +  ar  -ae  =  0  ’ 


r  — —  +  — -  +  2t  =0 
r  Sr  36  ^Tr6  * 


Also 


(2.6) 


where 


is  the  hydrostatic  part  of  the  stress  tensor,  and  s^  is  the  deviatoric 
part.  Using  equation  (2.6),  the  equilibrium  equations  reduce  to 


.  3s  n  3  s  ^  s  -  So 

3p  _  +  1  r9  +  _I _ _£ 

3r  "  3r  r  36  r 


l£  = 

36 


r6 

3r 


3s 

— —  +  2s 

36  r8* 


which,  upon  noting  that 


become 
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8r 


1  3sre 
r  80 


(2.7) 


^  =  -  !!n  +  2s  +  r 
80  80  br0  r  8r  * 


In  terms  of  the  strain  rates,  given  by  equations  (2.5),  the  invariant  I 
becomes 


p  p  1/2  -I  o  p  1/2 

I  =  [4drZ  +  4dr02]  =  -j  [4g2  +  (g1)2]  .  (2.8) 

Substitution  of  the  relation  (2.8)  and  the  strain  rates  (2.5)  in  the 
constitutive  equations  (1.4)  gives  for  the  components  of  the  stress 
deviators 


sr  =  -  2kg/[4g2  +  (g1)2] 


1/2 


sr0  =  ^  +  kg'/[4g2  +  (g 1  )2] 


1/2 


(2.9) 


The  hydrostatic  pressure  p  can  be  eliminated  from  equations  (2.7),  and 
the  final  result  is 


2 


02S. 


8r80 


2  9sr  ,  1  5%0 
r  80  r  902 


82s 


8r 


r0 

T 


8s 


-  3 


r0 


8r 


=  0. 


(2.10) 


Upon  substitution  of  the  stress  deviators  (2.9)  in  equation  (2.10),  an 
ordinary  differential  equation  is  obtained  for  g(0);  the  differential 


equation  is 


Mg'" 

r 


+  4g 1 ]  +  k[{ 


[4g2  +  (g')2]1/2 


I  I 


4{ 


[4g2  + 


(g 


)2]1/2 


-}  ]  =  0. 
(2.11) 


Clearly  g  =  const,  satisfies  equation  (2.11),  consequently  this 
solution  cannot  depend  on  y  or  k  and  the  corresponding  flow  is  axially 
symmetric.  Alternatively  in  order  that  the  differential  equation  (2.11) 
be  satisfied,  it  is  therefore  necessary  that  g ( 6 )  satisfies  both 


g"  '  +  4g'  =0  , 


(2.12a) 


and 


{ 


[4g2  + 


3  21/2> 

(g1)2]  /2 


i  i 


-  4{ 


[4g2  + 


.2 


(g'j*?7* 


•}  =  0. 


(2.12b) 


It  is  noted  that  g ( 0 )  is  governed  by  two  ordinary  differential  equations, 
whereas  for  the  velocity  field  obtained  by  Hill  (30)  for  the  plane  rigid 
perfectly  plastic  solution,  there  is  no  r  dependence,  and  g ( 0 )  is  defined 
by  one  equation. 

Further  constraints  on  g ( 0 )  are  obtained  from  the  boundary 
conditions  at  6  =  0,  6  -  a,  and  from  a  specification  of  the  volume  flow 
per  unit  length. 

From  symmetry 


g'(0)-0.  (2.13) 

Further,  since  the  Bingham  solid  may  be  regarded  as  a  viscous  medium  with 
a  variable  coefficient  of  viscosity,  it  appears  reasonable  to  impose  the 
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no-slip  boundary  condition,  usually  assumed  for  viscous  flow,  that  is 

g(a)  =  0. 


If  the  volume  flow  per  unit  length  is  Q,  then 


Q  =  2 


a 


g(s)  ds 


The  solution  to  the  differential  equation  (2.12a)  is  readily 
found  to  be 


g ( 0 )  =  A  +  B  cos  26  +  C  sin  20, 

where  A,  B,  and  C  are  arbitrary  constants  of  integration.  The  symmetry 
condition  (2.13)  requires  that 


consequently 


g(0)  s  A  +  B  cos  20.  (2.14) 

Substitution  of  equation  (2.14)  in  the  differential  equation  (2.12b) 
yields 

4B  sin28[-3A2B2  cos22e  -  2(A3b  +  4AB3)cos2e  -  (A4  +  2B4)]  =  0.  (2.15) 
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Equation  (2.15)  possesses  the  real  and  trivial  solution 

B  =  0. 


Further,  imposition  of  the  boundary  condition  at  0  =  a  implies  that 

A  =  0. 

Thus  the  assumption  of  radial  flow,  subject  to  the  no-slip  boundary 
condition  at  6  =  a,  has  yielded  the  trivial  solution. 

A1 ternatively ,  one  might  attempt  to  employ  the  solution 

g(e)  =  d, 


in  the  region 


0|£!S;<|ai  , 


where  D  is  a  constant,  and  the  solution 


g ( 0 )  =  A  +  B  cos  26  +  C  sgne  sin  26 


in  the  region 


3  <  0  <  a 


■* 


' 
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subject  to  the  constraints  noted,  plus  the  requirements  that  g ( 6 )  and 
g'(6)  be  continuous  at  0  =  ±  3. 

Continuity  of  g 1 ( 0 )  at  0  =  3,  therefore,  requires  that 

g ' (3)  =  o. 


This  condition  requires  that 


C  =  Bsgn  3tan  23.  (2.16) 

In  order  to  retain  continuity  in  the  velocity,  and  hence  in  g ( 0 ) ,  at 
0  =  3,  it  is  required  that 

D  =  A  +  B  cos23  +  Csgn  3sin  23, 

or  using  relation  (2.16) 


D 


A  +  — - — 

M  cos23  * 


consequently 


g ( 0 )  =  A  +  B  [cos20  +  tan23sin20].  (2.17) 

The  expression,  which  is  obtained  by  substitution  of  equation  (2.17) 
in  the  differential  equation  (2.12b),  is  very  complex  and  is  not  reproduced 
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here.  It  can  be  shown  that  the  non-trivial  solution  to  this  expression 
results  in 


B  =  0. 


This  case,  however,  has  been  considered  before  and  gives  rise  to  the 
trivial  solution,  with  zero  radial  velocity.  Hence  a  solution  of  the  as¬ 
sumed  form  does  not  exist. 


2.2  Axially  Symmetric  Visco-Plastic  Flow 

This  section  deals  with  the  axially  symmetric  converging  flow  of 
an  incompressible  Bingham  solid  in  a  circular  cone.  The  z-axis  is  taken 
along  the  axis  of  symmetry  and  spherical  polar  coordinates,  ( r ,4> ,0 )  •  are 
used  in  considering  the  flow.  The  flow  is  assumed  to  be  radial,  steady, 
and  quasi-static.  It  is  further  assumed  that  the  cone  is  long  enough  to 
neglect  any  end  effects.  The  cone  may  be  described  as  the  part  bounded 
by  <J>  =  a  and  <J>  =  -  a,  where  a  is  the  semi-angle  of  the  cone,  see  FIG.  2.2. 

Denoting  the  outward  radial  velocity  by  v^,  then  by  the  velocity 
strain  tensor  (2.1),  the  only  non-vanishing  strain  rates  are  given  by 
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FIGURE  2.2 

Visco-Plastic  Flow  in  a  Circular  Cone 


The  continuity  equation, 


9v 


r 


9r 


+ 


requires  that 


v 


r 


(2.19) 


Substitution  of  this  velocity  field  (2.19)  in  the  strain  rate  relationships 
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(2.18)  gives 


,  where 


d_ 

”  d<J>  * 


(2.20) 


The  non-vanishing  components  of  the  stress  tensor  are  ar,  aQ,  and  x^. 
Neglecting  inertia  terms  and  body  forces,  the  equations  of  equilibrium  are 


do  1  9l  -I 

7—^-  + - rr*’  +  —  (2a  -  a,-aQ  +  T  cot<(>)  =  0, 

9r  r  94>  r  r  <j)  6  r<j>  y 


[(a<0  -  °e)cot*  +  3Tr*]  =  °' 


Substitution  of  relation  (2.6)  in  the  equations  (2.21)  yields 


(2.21) 


~ -  +  - — +  —  (2s  “  s .  -  sA  +  s  .  cot  4>) , 

9r  9r  r  9<J>  r  '  r  <f>  6  r4>  * 


H  =  r  -fei  +  ir+  <s*  -  se)  cot*  +  3sr<J>  > 


(2.22) 


’6 


=  “  (s, 


V 


and  upon  noting  that 
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equation  (2.22)  becomes 


= 

9r 


+ 


1_ 

r 


sr(f)COt^* 


9<J> 


9s  , 

— £  +  3S 

94)  sr4>* 


(2.23) 


In  terms  of  the  strain  rates,  as  given  by  equations  (2.20),  the  invariant 
I  becomes 

I  =  [2dr2  +  2^*  +  2d02  +  4d^2]V2  =  -^j[12g2  +  (g')2]V2  .  (2.24) 

Substitution  of  the  relation  (2.24)  and  the  strain  rates  (2.20)  in  the 
constitutive  equations  (1.4)  gives  for  the  components  of  the  stress  deviators 

s  =  -  M  .  4kg/[12g2  +  (g')2]V2  , 
r 

1/2 

s4>  =  s0  =  3  +  2k3/[1292  +  (g')2]  .  (2.25) 

.  0  0  1/2 
S  =  3~+  kg,/[12g  +  (g,)2]  • 

The  hydrostatic  pressure  p  may  be  eliminated  from  equations 
(2.23),  and  the  resulting  partial  differential  equation  in  terms  of  the 
stress  deviators  is 


3%.  +  3^, 

9r94>  r  3cj> 


.  1 

r  ^2 


92s 


r4) 


9r 


r 


+ 
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+  fl?^cot*]  -  4^=  0.  (2.26) 

Substitution  of  the  stress  deviator  components  (2.25)  in  the  equation 
(2.26)  results  in  the  following  ordinary  differential  equation  for  gU): 


[g  1  '  1  +  (g '  co tcj) ) 1  +  6g '  ]  + 
r 


+  k[{ 


[12g2  +  (g')]1/2 


.}"  +  {  g'cot  0  } 


[I2g2  +  (g')2]1/2 


-  12{ 


[12g2  Mg')2]172 


■}']  -  0, 


Again  g  =  const,  satisfies  this  equation  and  the  flow  corresponding  to  this 
solution  is  spherically  symmetric.  Alternatively  it  follows  that 

g1  "  +  (g'cotcj)) '  +  6g '  =0, 


which  upon  integration  gives 


g '  1  +  g  1  cotcf)  +  6g  =  const  , 


(2.27a) 


and 

{ 


I  I 

+{ 


g 1  cot  <j) 


}  -  1 2{- 


[12g2  +  (g 1  )2]1/2  D2g2  +  (g')2]1/2  [12g2  +  (g')2]l/2 


r)  =  0. 


(2.27b) 


Again,  as  is  the  case  for  the  plane  flow  problem,  g(<j>)must  obey  two 
ordinary  differential  equations. 

The  condition  that  the  flow  be  symmetric  with  respect  to  the 
cj)  =  0  axis  gives 


g ' (o)  =  o. 


(2.28) 


29v f q  norJ6'i|)9in  r  noqu  rbMw 
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As  was  pointed  out  in  the  previous  section,  a  logical  boundary  condition 
to  assume  at  cf>  =  a  is 


g(a)  =  0  . 


(2.29) 


Furthermore,  the  volume  flow  condition  is 


Q  =  2tt 


a 

f 

o 


g(<j))  s i ncf)  dcf>. 


(2.30) 


The  problem  of  finding  the  unknown  function  g(<J>)  then  reduces  to  solving 
the  ordinary  differential  equation  (2.27a),  subject  to  the  differential 
equation  (2.27b),  and  boundary  conditions  (2.28)  and  (2.29),  and  volume 
flow  condition  (2.30).  The  solution  to  the  differential  equation  (2.27a) 
is  made  up  of  two  parts,  the  particular  solution 


gp(4>)  =  constant, 

and  the  homogeneous  solution,  g^tcj)),  which  satisfies 

gH ' '  +  gH'  cot<j>  +  6gH  =  0.  (2.31 ) 


Under  the  substitution 


x  =  COScf), 
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equation  (2.31)  becomes 

(1  -  x2)gH"  -  2xgH‘  +  2(2  +  l)gH  =  0,  (2.32) 

where  1  5  ^  ,  and  is  recognized  as  Legendre's  ordinary  differential 
equation.  The  solution  to  equation  (2.32)  may  be  written  in  the  form 

gH(x)  s  A  P2(x)  +  B  Q 2 ( x ) , 

where  P2(x)  and  Q2 ( x )  are  the  Legendre's  functions  of  degree  2,  of  the 
first  and  second  kinds,  respectively ,  and  A  and  B  are  arbitrary  constants 
of  integration.  Returning  to  the  <J>  variable,  this  solution  reduces  to 

gH(4>)  =  A  P2(cos<t>)  +  B  Q2(cos<|>)  , 


where 


P2(cos<J>)  =  j(3cos24)  -  1 ) , 


and 


Q2(cos<J>)  =  -l(3cos2cf>  -  1)  log  ]  *  .  1  cosep  . 

The  general  solution  to  the  differential  equation  (2.27a), 

g(<J>)  =  9p(4>)  +  9H(0)  » 


,(x)<.0  8  ^  (x)»q  a  =  (x)4e 
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then  reduces  to 

g(4>)  =  A(3cos2(j)  -  1)  +  B[ ( 3cos24)  -  1 ) T og  ]  *  -  6cos4>]  +  C, 

where  A,  B,  and  C  are  to  be  determined  such  that  the  boundary  conditions 
(2.28)  and  (2.29),  and  the  differential  equation  (2.27b)  are  satisfied, 
and  further  that  g(<|>)  be  finite. 

The  latter  condition  requires  that  B  =  0.  The  symmetry  condition 
is  satisfied  automatically.  There  remains  the  problem  to  find  A  and  C 
such  that  the  differential  equation  (2.27b)  is  satisfied  along  with  the 
boundary  condition  (2.29).  Substitution  of 

g(4>)  =  A(3cos2<t>  -  1 )  +  C 

in  the  differential  equation  (2.27b)  gives 

64Asin<t>cos<)>[(68A3C  -  22A4)cos6<f>  +  (348A2C2  +  7A4  +  U1A3C)cos4* 

(2.33) 

+  (A3C  +  318AC3  +  1  53A2C2)cos2<J>  +  (108C4  +  18A2C2  +  111AC3)]  =  0. 

The  real  and  non-trivial  solution  to  equation  (2.33)  gives 

A  =  0. 


The  solution  common  to  both  differential  equations  therefore  reduces  to 
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g(4>)  =  c, 

and  the  boundary  condition  (2.29)  requires  that  C  be  zero.  Consequently 
the  assumption  of  radial  flow  subject  to  the  no-slip  boundary  condition 
at  <J>  =  a  has  yielded  the  trivial  solution  in  both  the  plane  and  axially 
symmetric  problems. 

It  can  be  shown  that  an  attempt  to  treat  the  flow  as  composed 
of  two  regions,  as  was  discussed  in  the  plane  flow  case,  again  leads  to 
the  trivial  solution. 

The  solution  of  both  visco-plastic  problems  clearly  involves 
curved  streamlines,  and  to  obtain  a  solution  it  is  probable  that  the  inlet 
(or  exit)  velocity  profile  must  be  known.  A  numerical  approach,  using  a 
relaxation  procedure,  was  attempted  but  was  abandoned  since  the  solution 
rapidly  diverged. 
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CHAPTER  III 


PSEUDO-PLASTIC  FLOW 

In  CHAPTER  II  the  plane  radial  flow  of  an  incompressible  Bingham 
solid  in  a  converging  channel,  and  the  axially  symmetric  converging  flow 
in  a  circular  cone  were  introduced.  It  is  recalled  that  g ( 0 )  in  both 
cases  was  constrained  by  two  ordinary  differential  equations.  It  was  shown 
that  under  the  no-slip  boundary  condition  only  the  trivial  solutions,  with 
zero  radial  velocities,  were  obtained.  It  was  then  suggested  that  a  more 
general  formulation  assuming  non-radial  flow  is  intractable. 

An  alternative  model,  the  pseudo-plastic  material,  is  now  con¬ 
sidered.  This  model  exhibits  rate  effects  and  under  the  assumption  of 
radial  flow  removes  the  r-dependence  in  the  governing  equations  and  yields 
one  ordinary  differential  equation  for  the  required  function.  The  absence 
of  the  r-dependence  is  similar  to  that  obtained  by  Nadai  (29)  and  Shield  (31) 
for  the  rigid  perfectly  plastic  materials. 

It  is  shown  in  the  body  of  this  chapter  that  the  boundary  value 
problem  for  the  required  function  can  be  reformulated  as  an  initial  value 
problem  for  a  third  order  non-linear  ordinary  differential  equation  in  g (0 ) . 
This  problem  is  readily  solved  numerically  by  employing  Gill's  variation  of 
the  Runge-Kutta  fourth  order  method. 

3. 1  Plane  Pseudo-Plastic  Flow 

The  plane,  radial,  steady,  quasi-static  flow  of  an  incompressible 
pseudo-plastic  material  in  a  converging  channel  is  considered.  FIGURE  3.1 
shows  some  necessary  notation.  As  previously  discussed,  the  channel  is  as- 
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sumed  to  be  sufficiently  long,  so  that  the  velocity  field  and  the  state 
of  stress  are  independent  of  the  conditions  occurring  at  the  entry  and 
exit  sections. 
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Pseudo-Plastic  Flow  in  a  Converging  Channel 


The  constitutive  equation  for  the  pseudo-plastic  material  is 


s 


ij 


n-1 


9 
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which  upon  applying  equations  (2.5), 


dr  = 


g(0) 

2 


d  = . d  = sill 
e  r  r2 


(2.5) 


yields 


d  =  all6-I 

re  TT  ’ 


3-n 

9  2 

s  =  _  c  =  _  2 _ k 

r  9  r2n 

1-n 


O  ^ 

 2  u 


n-1 

{g[4g2  +  (g1)2]  2  }  , 


n-1 


sr6  =  "2F"  {9'  [4g2  +  (9')2;l  2  K 


Let 


g ( 9 )  =  Am(e)  ,  (3.1) 

where  m(6)  is  a  non-dimensional  function  of  6,  and  A  is  a  constant  to  be 
determined  by  the  volume  flow  per  unit  length  condition 


m(e)  de }  . 


Using  equation  (3.1),  the  stress  deviator  components  become 
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3-n 


s 


r 


1-n 


(3.2) 


where 


n-1 


6  =  [4m2  +  (m ' ) 2]  2 


The  differential  equation  governing  m(0)  is  obtained  by  substitution  of 
equations  (3.2)  in  equation  (2.10),  and  is 


(m 1 G) ' '  +  4n(l-n)m'G  +  4(2n-l)(mG)'  =  0  . 


(3.3) 


It  is  significant  that  there  is  no  r-dependence  in  equation  (3.3),  unlike 
the  governing  differential  equations  for  the  problems  considering  the 
Bingham  solid.  Upon  replacement  of  the  stress  deviator  components  in 
equations  (2.7)  by  the  relations  given  by  equations  (3.2),  the  expression 
for  the  hydrostatic  pressure  p  is  found  to  be 


1-n 


P  -  [f(l-n)  mG  -  ^  (m'G) 1  ]  +  C, 


where  the  constant  C  depends  on  the  prescribed  pressure  at  some  point  in 


the  flow  field. 

If  n  =  0,  equation  (3.3)  reduces  to  the  equation  for  Nadai's 


norJi.  rUadua  ^d  b9nUJdo  zt  (e)m  nh;i9 
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problem*,  and  n  =  1  yields  the  governing  differential  equation  for  the 
creeping  flow  of  a  Newtonian  viscous  liquid. 

The  equation  (3.3)  appears  to  be  in  general  analytically  in¬ 
tractable,  but  the  case  n  =  j  can  be  integrated  and  becomes 


*Substi tuting 


reduces  equation  (3.3)  to 


n  =  0 


4[y- . -m~  -] 1  •  0  , 

/4m2  +  (m')2 


and  with  the  change  in  the  dependent  variable 


m 

2m 


this  becomes 


h  1  1  1  1 

[-===7-]  -  2  [— =r-]  =  0  . 

A  +  A  +  h2 


With  the  substitution 


h  =  -  tan  2\p , 

the  differential  equation  for  becomes 

ip '  -  c  sec  2\p  -  1 , 

where  c  is  a  constant  of  integration,  and  the  corresponding  differential 
equation  for  m  is 


ST  =  *  tan  2*  ' 

It  is  readily  verified  that  the  two  ordinary  differential  equations  for  \p 
and  m  are  the  same  as  those  obtained  by  Nadai  (29)  and  Hill  (30)  respectively, 
thus  illustrating  the  third  axiom  of  Rheology. 


o)  {ZA  no  Jsuo.  2 do u  > sn 


i r  dt' 1  Mj 


- 
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m'G  =  Bcos6  +  Csin6,  (3.4) 

where  B  and  C  are  constants  of  integration.  Since 

m ' (0)  =  0, 


equation  (3.4)  reduces  to 


m'G  =  CsinG, 


and  the  shear  stress  becomes 

1  1 

A  *2 

sr0  ■  C  sine  . 

This  expression  for  srQ  is  used  to  estimate  the  accuracy  of  the  numerical 
procedure  used  to  integrate  equation  (3.3)  for  arbitrary  n. 

The  boundary  value  problem  now  involves  the  determination  of  m(0) 
e  C3>  satisfying  equation  (3.3)  in  the  region  bounded  by  6  -  ±  a,  and 
subject  to  the  symmetry  condition 

m'(0)  =  0  , 

and  the  assumed  no-slip  boundary  condition 


m(a)  =  0  , 
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and 


m(0)  =  1, 


where  A  is  obtained  from  the  volume  flow  per  unit  length  condition 


Q  =  A{2 


a 


0 


m(e)  de}  . 


Considerable  effort  was  expended  on  unseccessful  attempts  to 
solve  this  boundary  value  problem  directly.  The  boundary  value  problem 
was  therefore  reformulated  as  the  following  initial  value  problem. 

Solve  equation  (3.3)  subject  to  the  initial  conditions: 


m ( 0 )  =  1, 


m 1  (0)  =  0 


from  symmetry,  and 


m'  1  (0)  =  3, 

where  3  is  determined  by  the  requirement  that 

m(a)  =  0, 


which  is  the  assumed  no-slip  boundary  condition.  It  is  known  that  such  a 
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problem  may  be  reduced  to  a  system  of  three  first  order  differential 
equations,  which  must  be  satisfied  subject  to  the  prescribed  initial 
conditions.  Reducing  equation  (3.3)  to  a  system  of  first  order  ordinary 
differential  equations  gives 


y*i  (e)  =  y2(e)  , 


y2  (e)  =  y3(e)  » 


y3  (e)  =  F[y i (e) ,  y2(e) ,  y3(e) ]  , 


where 


y-j (e)  =  m(e)  , 


y  2  ( e )  s  m '  ( e )  , 
y 3 ( e )  s  m' 1 ( e)  , 


and 

F[y,  (0)  ,y?(e)  ,y,(6)] - ] - j  {(4n2-12n+4)m'M 

-  [(4n2-6n+2)m  +  (n-l)m' ']M'  -  (n-1 )(n-3)m' (M‘ )2/4M 
-  (n-1  )m ' [4(m ' )2  +  4mm"  +  (m")2]}  , 


9teriw 
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where 


M  =  4m2  +  (m')2  . 


The  initial  values  become 

y,(o)  =  l, 
y2(o)  =  o, 
y3(o)  =  e. 

This  formulation  is  readily  amendable  to  numerical  treatment  using  Gill's 
variation  of  the  Runge-Kutta  fourth  order  method. 

The  method  of  solution  is  similar  to  that  used  for  the  axially 
symmetric  converging  flow  in  a  circular  cone,  and  the  details  are  included 
in  the  discussion  of  that  problem. 

The  stress  and  flow  fields  obtained  for  the  plane  flow  case  are 
applicable  to  the  compression  of  a  pseudo-plastic  material  between  two 
inclined  plates.  Letting  U  be  the  velocity  of  the  plates  normal  to  their 
lengths,  see  FIG.  3.2,  then  the  velocity  field  is  obtained  by  superimposing 
a  velocity  U/sina  to  the  diverging  flow  field  obtained  above.  The  corres¬ 
ponding  stress  field  is  that  obtained  for  diverging  flow.  This  solution 
will  be  valid  if  the  plates  are  large,  and  the  flow  is  considered  far 
enough  away  from  the  ends,  so  that  end  effects  do  not  influence  the  solution. 


2'  nort'jroe  o  borlJofn  srIT 
fu;  a  .  o  ^  ?*'  JVf  3 
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FIGURE  3c 2 

Compression  of  a  Pseudo-Plastic  Material 
Between  Inclined  Plates 

3.2  Axially  Symmetric  Pseudo-Plastic  Flow 

The  problem,  now  considered,  is  that  of  the  axially  symmetric, 
radial,  steady,  quasi-stati c ,  converging  flow  of  an  incompressible  pseudo¬ 
plastic  material  in  a  circular  cone.  Some  pertinent  notation  is  shown  in 
FIG.  3.2.  As  in  the  previously  discussed  flow  problems,  the  cone  is  as¬ 
sumed  to  be  long  enough  so  that  the  end  effects  may  be  neglected.  The 
non-vanishing  strain  rates,  for  radial  flow,  are 


t 
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FIGURE  3.3 

Pseudo-Plastic  Flow  in  a  Circular  Cone 


d 


rep 


A-  ,  where 
2r 


d_ 

d4> 


(2.20) 


Substitution  of  equations  (2.20)  in  the  constitutive  equations  (1.6)  results 
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in 


5-n 
9  2 

s  =  -  2 _ U. 

r  r3n 
3-n 

9  2 

S  =  S  =  - * - '£• 

4>  e  r3n 

1-n 


n-1 

2 

{g[12g2  +  (9')2]  }, 


n-1 

{g[12g2  +  (g')2]  2  }  , 


n-1 


sr<j>  =  ^^{g'E^g2  +  (g'}2]  2  }. 


(3.5) 


Substituting  in  equations  (3.5) 


g(cf))  =  Am(<J))  , 


where  m(4>)  is  a  non-dimensional  function  of  c|> : 


5-n 

.  2  2  uAn  - 

sr  3n  mG  ’ 

r 

3-n 

o  2  «n 

s4>  ‘  se  *  r3n  mG  ’ 
1-n 

o  2  An 

c  =  2  nA  m  i  r 

sr<(>  3n  m  G  • 

r 


(3.6) 


where 


n-1 

2  2  2 

G  =  [12rrr  +  (m‘ )  ] 


and  A  is  a  constant  to  be  determined  from  the  volume  flow  condition 


-n 

f-('m)  +  iSf]  0 
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rOL 


A{2tt 


0 


m( 4>)  sin<j>  dcf>}  . 


Substitution  of  equations  (3.6)  in  equation  (2.26)  gives 


(m'G) '  '  +  (m'Gcotcf)) '  +  9n(l-n) (m'G)  +  6(3n-2)  (mG) '  =  0  .  (3.7) 


Again  the  governing  differential  equation  has  no  r-dependence. 

After  substitution  of  equations  (3.6),  the  integration  of 
equations  (2.23)  results  in 

1-n 

2  n 

p  =  - — 3n^A  [pj-  (l-n)mG  -  ^  {(m'G)'  +  m'Gcot<}>}]  +  C  ,  (3.8) 

where  the  constant  C  depends  on  the  prescribed  pressure  at  some  point  in 
the  flow  field. 

As  a  check,  it  is  noted  that  the  substitution  n  =  0  in  equation 
(3.7)  yields  the  equation  obtained  by  Shield  (31)  for  the  perfectly 
plastic  von  Mises  solid*,  and  the  substitution  n  =  1  reduces  equation  (3.7) 
to  that  for  the  creeping  flow  of  a  Newtonian  viscous  liquid. 


*Substi tuting  n  =  0  in  equation  (3.7)  gives 

m  1  I  I  m  I  I  m  I 

-  ■  ■  ]  +  [  m  -  cot<|)]  -  12[._ - J1 ]  =  0,(3. 9) 

/l2m2  +  (m')2  /12m2  +  (m')2  /l2m2  +  (m')2 

and  with  the  change  in  the  dependent  variable 


r 


46 


equation  (3.9)  becomes 


?  1/2 

+  TCOt<|>  ±  2/3  (1-t  ) 


C, 


where  c  is  an  arbitrary  constant  of  integration.  Integrating  the  above 
expression  for  m  in  terms  of  t  gives 

2  -1/2 

m  =  B  exp  {±  2/3  /  t(1-t  )  d4>) . 

The  differential  equation  for  t,  and  the  expression  for  m  in  terms  of  x 
are  consistent  with  those  obtained  by  Shield  (31),  again  illustrating  the 
third  axiom  of  Rheology. 


Returning  to  a  consideration  of  equation  (3.7);  the  boundary 
conditions  to  which  the  solution  to  this  differential  equation  shall  be 
subjected  to,  are 


m(0)  =  1, 


m '  ( 0 )  =  0 


from  symmetry ,  and 


m(a)  =  0, 

that  is  the  no-slip  boundary  condition  at  the  walls. 

If  n  =  2/3,  the  differential  equation  reduces  to 

(m'G)  +  (m'Gcot4>)  +  2m'G  =  0  , 


which,  under  the  substitution 


,0  ?  (jd  )  nri 
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x  =  cos 


becomes 


p  1 1  i  ]2 

(1-x  ; (m'G)  -  2x(m'G)  +  [1(1+1) - *-]  m'G  =  0  ,  (3.10) 

1-x^ 


which  is  Legendre's  equation  of  degree  1  and  order  1.  The  solution  to 
equation  (3.10)  is 


m'G  =  BP1 1 (x)  +  CQ1 1 (x)  , 


where  pJ  and  qJ  are  the  associated  Legendre's  functions  of  degree  1  and 
order  1,  of  the  first  and  second  kinds,  respectively,  and  B  and  C  are 
arbitrary  constants  of  integration.  Analytically 


m 

p>)  =  (l-x2)2  ^*1  . 
n  dxm 


The  general  solution  to  the  differential  equation  (3.10)  therefore  becomes 


m'G  =  BP-^tcosef))  +  CQ-j1  (coscf>)  . 

The  symmetry  condition  requires  that  m'G  be  zero  at  4>  =  0;  since  Q^(0) 
is  not  bounded  at  <j>  =  0,  it  is  necessary  that 


1/ 


.r  'teb'to  b n*>  f  stngsb  nottsups  a  yibnsgsJ  af  rtorriw 
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C  =  0, 


consequently 


m'G  =  Bsincj)  . 


(3.11) 


Applying  equation  (3.11),  the  shear  stress  becomes 


1  2 


(3.12) 


It  was  not  possible  to  obtain  a  closed  form  solution  for  other 
values  of  n  in  the  range 


0  <  n  <  1  , 


and  equation  (3.12)  was  used  to  estimate  the  accuracy  of  the  numerical 
procedure  used.  The  boundary  value  problem  can  be  more  readily  solved  by 
transforming  it  into  an  initial  value  problem. 

This  initial  value  problem  requires  the  solution  of  equation 
(3.7)  subjected  to  the  initial  values 


m(0)  =  1 , 


m'(0)  =  0 


from  symmetry,  and 


Jri9Up92flOD 


.  r  .  ■  .s 
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m' 1 (0)  =  3. 

The  parameter  3  is  determined  by  the  no-slip  boundary  condition 

m(a)  =  0. 

Letting 

y-j  (4>)  =  m(4>)  , 
y2(<J>)  =  rn '  (4>)  , 
y3(<l>)  =  rn'  '(<(>)  , 

then  equation  (3.7)  is  equivalent  to  the  system 

i 

y-j  (4>)  =  y2U)  , 

t 

y2  (<J>)  =  y3(4>)  , 

y3  U)  =  F[<|> ,y-j (<J>)  •  y2(4>) »  y3U)L 

where 

F[<t>»y1  4)  .y2(<t>)  .y3(<}>) ] - ] - J  {[(9n2-27n+12)m' 

M  +  (n-1 ) (m 1 ) 


r,  -  r:>!  )3  y.11  Muod  qHz-o-  a,'.'  \  t-  m  -f  8  •>  '  I  9r  r 


,  U)  'm  =  (o>)r\6 
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-  cotym"  +  csc2<J)m']M  -  [(9n2-15n+6)m  +(n-l)m" 
+  (^)cotcf>m']  M'  -  (n-l)(n-3)m'(M')2/4M 
-  (n-1  )m ' [1 2(m '  )2  +  12mm1'  +  (mM)2]}, 
and 


M  =  12m2  +  (m')2  . 

Gill's  variation  of  the  Runge-Kutta  fourth  order  method  is  used 
to  solve  the  initial  value  problem.  The  precise  procedure  utilized  is 
considered  briefly  in  the  next  section. 

After  this  work  had  been  nearly  completed,  the  author  was  shown 
a  paper  by  Tanner  (32),  in  which  equation  (3.7)  was  presented.  The  nu¬ 
merical  method  of  solution  is  not  indicated  in  this  paper,  and  it  does  not 
include  the  partial  analytical  solution  for  the  special  case  with  n  =  2/3. 

3.3  Integration  Procedure 

It  has  been  shown  that  both  the  plane  and  the  axially  symmetric 
boundary  value  problems  are  reducable  to  equivalent  initial  value  problems 
involving  the  solution  of  a  system  of  first  order  ordinary  differential 
equations,  subjected  to  prescribed  initial  conditions. 

The  initial  conditions  to  be  Imposed  on  m( 0)  and  m'(e)  are  clear; 
however,  the  initial  value  of  mM(e)  is  not  readily  prescribed.  Rather, 
a  trial  and  error  method  is  adopted;  various  values  of  mM(0)  are  assumed 
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until  the  desired  boundary  condition  of  m(a)  =  0  is  obtained. 

Gill's  variation  of  the  Runge-Kutta  fourth  order  method  is 
employed  as  an  integration  procedure.  An  explanation  of  this  method, 
and  a  brief  discussion  on  its  accuracy,  is  contained  in  the  Appendix. 

The  Appendix  also  contains  a  copy  of  the  DRKGS  subroutine.  This  is  the 
IBM  System/360  Scientific  Subroutine  Package  for  Gill's  variation  of  the 
fourth  order  Runge-Kutta  method. 

It  has  already  been  noted  that  possession  of  a  closed  form 
solution  for  s^,  with  n  =  2/3,  for  the  axially  symmetric  problem,  and 
for  srQ,  with  n  =  1/2,  for  the  plane  problem  permits  a  direct  comparison 
of  the  numerical  solution  to  the  closed  form  solution,  and  hence  an  assess¬ 
ment  of  the  accuracy  can  be  made  for  other  values  of  n. 

The  values  of  n  considered  are 

0  <  n  <  1 . 

The  value  n  =  0  corresponds  to  the  perfectly  plastic  von  Mises  solid, 
and  n  =  1  to  the  Newtonian  viscous  liquid;  hence,  it  is  possible  to  study 
the  transition  from  the  perfectly  plastic  von  Mises  solid  to  the  Newtonian 
viscous  liquid. 

The  results  of  the  analysis  performed  in  this  chapter  are  presented 
in  FIGS.  6.II-1  to  6. 11-13  for  the  plane  problem,  and  in  FIGS.  6.III-1  to 
6.III-15  for  the  axially  symmetric  case,  and  a  brief  discussion  of  these 
results  is  given  in  CHAPTER  VI. 
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3.4  On  the  Inclusion  of  Inertia  Effects 

The  solutions  presented  in  this  chapter  for  the  flow  of  pseudo¬ 
plastic  materials  are  based  on  the  assumption  that  the  non-linear  term  in 
the  equilibrium  equation  can  be  neglected. 

The  inclusion  of  the  inertia  term  is  now  considered  briefly. 
Instead  of  the  equilibrium  equations  (2.7),  the  equations  of  motion  are 
used,  which  for  steady  flow  become 


3r 


+  pvr 


1 

r 


. 

30 


+  2s 


re 


+ 


(3.13) 


where  p  is  the  density.  Substitution  of  the  stress  deviators  (3.2)  in 
equation  (3.13) ,  gives 

l~n  2  2 

^  9  2  «n  pA  m 

3^r  =  — 2n+'l  [^n(mG)  -  4mG  +  (m1  G) 1  ]  +  — — >  (3.14) 

1-n 

|£  =  2  ^fl  [2(mG) '  +2m'G  -  2nm'G].  (3.15) 


From  equations  (3.14)  and  (3.15)  it  follows  that 


1-n 

2  2  uAn  r4mG  -  4n(mG)  -  (m '  G) 1 )  -, 

p  '  „2n  L  2n  J 

r 


2  2 
pA  m^ 

2 

2r 


C, 


(3.16) 


where  the  constant  C  depends  on  the  prescribed  pressure  at  some  point  in 
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the  flow  field,  and  further, from  equation  (3.16),  that 

1-n 

.  2  2  ufln  r4(mG) '  -  4n(mG) '  -  (m'G)"-.  pA2mm ' 

36  Tn  L  2n  J  “  2 

r  r 

Equating  equations  (3.15)  and  (3.17)  gives 


1-n 

2  n  2 

—  ■$-  {(m'G)"  +  4n  ( 1  -n )  (m 1 G)  +  4(2n-l )  (mG) ' }]  +  pA  =  0. 

(3.18) 


Consequently,  for  n  f  1 


(m'G)"  +  4n  ( 1  -n )  (m'G)  +  4(2n-l )  (mG) '  =  0, 


(3.19) 


and 


P—P  ■  =  0.  (3.20) 

r 

Equation  (3.20),  and  the  no-slip  boundary  condition,  imply  that 

m  =  0  . 

Equation  (3.19)  is  the  same  as  that  obtained  if  the  inertia  term  is  neglected, 
see  equation  (3.3) . 

It  is  interesting  to  note  that  if 
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n  =  1 


» 


then  equation  (3.18)  implies  that 


m ' 1 '  +  4m 


+  pAmm1  =  0 

y 


since 


G|  =  1  ; 
n=l 

this  is  the  governing  equation  for  the  plane,  steady  flow  of  a  Newtonian 
viscous  liquid,  without  body  forces,  in  a  converging  channel,  which  was 
considered  in  detail  by  Rosenhead  (26). 

Similarly,  the  inclusion  of  the  inertia  term  in  the  axially 
symmetric  problem  shows  that  there  is  no  radial  flow  solution,  not  even 
for  n  =  1 . 

3.5  On  the  Neglect  of  Inertia  Effects 

It  was  shown  above  that,  under  the  prior  assumption  of  radial, 
steady  flow,  the  presence  of  the  non-linear  inertia  term  in  the  equation 
of  motion  yielded  the  trivial  solution  with  zero  radial  velocity,  both 
for  the  plane,  if  n  t  1,  and  the  axially  symmetric  pseudo-plastic  flow, 
at  least  if  the  material  is  assumed  incompressible,  and  without  body  forces. 

Solutions  are  obtained,  both  for  the  plane  and  the  axially 
symmetric  flow,  with  the  assumption  that  the  flow  is  quasi -stati c.  This 
assumption  is  valid  if  certain  conditions,  which  relate  the  inertia  force 
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to  the  dissipative  forces,  are  satisfied. 

A  comparison  of  the  relative  magnitudes  of  the  inertia  term 
with  those  of  the  terms  retained  in  the  equilibrium  equation,  leads  to 
inequalities  which  are  of  little  value  since  numerical  results  make 
physical  interpretation  of  these  conditions  difficult.  A  better  approach, 
due  to  Batchelor  (33),  who  considered  the  inertia  effects  for  the  steady, 
quasi-static  flow  of  a  Newtonian  viscous  liquid  in  slowly-varying  circular 
cones,  can  be  applied  to  the  pseudo-plastic  flow  problems  considered  in 
this  thesis.  This  method,  however,  requires  a  knowledge  of  the  solution 
for  the  flow  of  the  pseudo-plastic  material  in  a  circular  tube,  which  is 
considered  first.  This  solution  is  then  used  to  specify  a  condition,  which, 
if  satisfied,  justifies  neglecting  the  inertia  term  for  the  axially  symmetric 
pseudo-plastic  flow  in  slowly-varying  circular  cones. 

Let  polar  coordinates,  (r,6,z),  describe  the  pseudo-plastic  flow 
in  the  circular  tube  of  radius  R;  the  z-axis  is  taken  as  the  axis  of 
symmetry,  see  FIG.  3.4.  The  flow  is  assumed  to  be  steady  and  quasi -static; 
furthermore,  it  is  assumed  that  the  velocity  profile  is  independent  of  6 
and  z.  The  incompressible  pseudo-plastic  material  under  consideration  is 
the  same  as  that  used  previously  in  this  thesis. 

Let  the  velocity  in  the  z-di recti  on  be  denoted  by  w,  which  by  the 
prior  assumption  is  a  function  of  r  only.  Consequently,  the  only  non¬ 
vanishing  component  of  the  strain  rate  tensor  becomes,  from  equation  (2.1), 

■  _  1  dw 

rz  1  dr  * 


(3.21) 
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FIGURE  3.4 

Pseudo-Plastic  Flow  in  a  Circular  Tube 


and  hence  the  strain  rate  invariant  I  becomes 


I 


dw 

dr 


The  only  non-vanishing  stress  tensor  component  is  srz,  which  by  substitution 
of  equation  (3.21)  in  the  constitutive  equation  (1.6)  becomes 

1-n 


0  2  /  dWxn 

srz  =  2  ' 


(3.22) 


Neglecting  body  forces,  the  equilibrium  equations  for  steady,  quasi¬ 


static  flow  are 
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9£  = 

3r 


0  , 


3£  =  !!n+  tlL=  li_  (rs  \ 

9z  9r  r  r  9r  '  rz; 
Equations  (3.23)  and  (3.24)  imply  that 


(3.23) 


(3.24) 


(3.25) 


P  =  p(z). 

Substitution  of  equation  (3.22)  in  equation  (3.25),  with  n  +  0,  results 
in 


dz 


1  d 

^  F  dr  ^dr^ 


n 


» 


since  the  left  hand  side  is  a  function  of  z  alone,  and  the  right  hand 
side  of  r  alone,  it  follows  that 


= 

dz 


1-n 


1  d 


y  r 


Tr(~) 
dr  Lndr; 


n 


]  =  c, 


(3.26) 


where  C  is  a  constant.  The  constant  C  is  associated  with  the  pressure 
difference  between  two  distant  sections,  A  and  B,  see  FIG.  3.4,  which 
causes  the  flow.  Substitution  of 
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dp  _  Ap 

dz  L  ‘  L 


in  equation  (3.26)  gives 


n-1 

[r(^.)n]  =  _  .2..„2  AP  r 

dr  Lr^dr  J  yL 


(3.27) 


Integrating  equation  (3.27)  once,  results  in 

n-1 

/ dw \  _  2  Ap  /A  r\ 

'dr'  "  yL  'r  “  2'  ' 


(3.28) 


where  A  is  a  constant  of  integration.  From  symmetry, 


dw  I  =0 

dr  I  u  > 

ar  r=0 


it  follows  that 


A  ■  0  . 


Consequently,  from  equation  (3.28), 


n-3 


1  1 


dw  92n  /AP\n  n 
dr  "  c  W  r 


(3.29) 


Integrating  equation  (3.29)  gives,  for  n  f  -  1, 


n-3  1 


1+n 


W  -  22n  (f )"  [B  -  O  r  "  ] 


The  parameter  B  is  a  constant  of  integration,  which  is  found  from  the 
assumed  no-slip  boundary  condition  at  the  walls,  and  is 
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1+n 

p  _  n  n  n 

B  ■  TTF  R 


The  velocity  profile  in  terms  of  the  pressure  difference,  therefore,  is 


n-3 


1  1+n  1+n 


tl  _  92n  /  n  wAPxn  ro  n  ^  n 

w  "  2  (TT7T)  (Tir)  [R  ‘  r  ]  * 


From  the  volume  flow  condition, 


rR 


Q  = 


2Trr  wdr  , 


the  velocity  profile  may  be  written  in  terms  of  Q  as 


1+n 

•  ■  ^  o  i’  -  ® "  i  ■ 


(3.30) 


It  may  be  seen  that  with  the  substitution 


n  -  1 


equation  (3.30)  reduces  to  the  case  for  the  Newtonian  viscous  flow. 

The  solution  obtained  is  now  used  to  obtain  a  condition,  which, 
if  satisfied,  justifies  neglecting  inertia  effects  for  the  axially  symmetric 
pseudo-plastic  flow  in  circular  cones  of  small  semi-angles.  The  stream- 
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lines  for  slowly-varying  circular  cones,  however,  are  not  uniaxial,  that 
is  parallel  to  the  z-axis,  as  is  the  case  for  the  flow  in  a  tube.  Hence, 
in  addition  to  the  uniaxial  velocity  w,  there  is  a  radial  velocity  u  which 
is  of  the  order  aw,  where 


is  the  inclination  of  the  streamline  to  the  z-axis,  see  FIG.  3.5. 


FIGURE  3.5 

Velocity  Components  in  a  Circular  Cone  of  Small 

Semi -Angle  a 


w 
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If  the  flow  in  a  circular  cone,  of  small  semi-angle,  is  considered, 
equation  (3.30)  is  a  valid  approximation  to  the  pseudo-plastic  flow  in 
that  cone,  provided  that  the  neglected  inertia  terms  are  small  compared 
to  the  term  retained  in  the  equilibrium  equation  (3.25). 

The  representative  magnitude  of  each  of  the  neglected  inertia 

terms , 


A11  dw  .  3w 
pu  w  ,  and  pw  ^ 


is 


pa 

The  representative  magnitude  of  the  term  retained  is 

qn 

Hence  equation  (3.30)  is  consistent  with  the  neglect  of  the  inertia  term 
if 


Q2 

pa  « 

K 


Qn 

V"+1 


> 


where  a  is  taken  as  its  maximum  value,  that  is  the  semi-angle  of  the  circular 
cone.  The  condition  specifying  the  justification  for  neglecting  the  in¬ 
ertia  term  therefore  becomes 
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p02"n  . 

a  «  1  • 

yR 


(3.31) 


Upon  introduction  of  the  characteristic  velocity  U,  defined  by 


u  = 

ttR^ 


where  R  is  the  characteristic  length  of  the  tube,  the  inequality  (3.31) 
reduces  to 


a 


Pu2-y 

y 


«  l  . 


(3.32) 


It  is  noted  that  for  n  =  1,  which  corresponds  to  the  Newtonian  viscous 
flow,  equation  (3.32)  reduces  to  that  obtained  by  Batchelor  (33).  Further¬ 
more,  it  is  observed  that  the  quantity 

pU2~nRn 

y 

in  equation  (3.32)  is  consistent  with  Reynold's  number,  as  defined  for 
the  power  law  model. 

Since  neglecting  of  inertia  terms  for  the  axially  symmetric 
pseudo-plastic  flow  problems  is  justified,  if  condition  (3.32)  is  satis¬ 
fied,  the  results  should  be  applicable  to  technological  processes  such  as 
wire  drawing  and  hydraulic  extrusion. 
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CHAPTER  IV 


APPLICATION  TO  WIRE  DRAWING  AND  HYDRAULIC  EXTRUSION 


In  CHAPTER  III  the  velocity  and  stress  fields  were  obtained 
for  the  plane  pseudo-plastic  flow  in  a  converging  channel,  and  axially 
symmetric  converging  flow  in  a  circular  cone. 

The  purpose  of  the  present  chapter  is  to  study  the  application 
of  the  axially  symmetric  solution  to  wire  drawing  and  hydraulic  extrusion 
the  working  material  shall  be  approximated  by  the  pseudo-plastic  model 
with  n  close  to  zero. 

FIGURE  4.1  indicates  that  the  pseudo-plastic  model, 


s 


•  • 
ij 


» 


with  n  =  0.1  should  be  a  suitable  model  for  a  rate  dependent  material. 
This  model  is  useful  for  application  to  extrusion  processes;  however,  it 
does  not  possess  a  yield  limit.  This  makes  it  difficult  to  apply  to  wire 
drawing  problems.  Since  wire  drawing  involves  uniaxial  tension  in  the 
drawn  material,  and  possibly  in  the  undrawn  material,  if  there  is  back- 
pull,  it  is  difficult  to  consider  this  process  if  there  is  no  yield  limit 
It  is  usually  assumed  that  the  yield  limit  is  not  exceeded  in  the  drawn 
and  undrawn  material,  but  only  in  the  deforming  material  passing  through 
the  die. 

The  following  procedure  can  be  adopted  to  circumvent  this  diffi 
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FIGURE  4.1 


Flow  Curves  for  the  Power  Law  Model 
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culty  and  to  enable  the  pseudo-plastic  model  to  be  applied  to  the  wire 
drawing  process.  Consider  the  minimum  value  I  of  I  in  the  region  ABCD, 
see  FIG.  4.3,  for  the  flow  field  obtained  from  equation  (3.7);  this 
minimum  occurs  at  point  E,  see  FIG.  4.3.  The  flow  curve,  see  FIG.  4.1, 
gives  at  the  point  I  =  IQ  the  corresponding  value  for  J/y,  denoted  by  J/y. 
The  assumption  is  now  made  that  the  material  stays  rigid  until  J  reaches 
its  critical  value  J*  <  JQ  .  The  mean  drawing  stress  must  be  less  than 
the  tensile  yield  stress,  which,  corresponding  to  the  above  assumption, 
is 


Y  =  /3  J*  <  /3  JQ.  (4.1) 

The  assumed  flow  curve  for  the  material  in  this  drawing  process,  then, 
is  as  shown  in  FIG.  4.2. 

Consider  FIG.  4.3,  which  represents  a  cross-section  of  a  conical 
reducing  die.  It  is  clear  that  severe  approximations  are  necessary  in 
order  to  analyse  these  processes.  The  method  used  in  this  thesis  shall 
follow  that  of  Shield.  Shield  (31)  presented  an  approximation  which  assumes 
that  the  stress  and  velocity  fields  of  the  material  in  the  region  ABCD, 
see  FIG.  4.3,  can  be  represented  by  those  obtained  in  CHAPTER  III;  this 
approximation  is  reasonably  valid  for  small  semi-angles  and  large  reductions, 
and  neglects  the  inlet  and  exit  effects,  sections  AB  and  CD  respectively, 
see  FIG.  4.3.  These  and  other  approximations  will  become  evident  in  the 
body  of  this  chapter. 

Theoretical  results  are  presented  for  both  wire  drawing  with 
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FIGURE  4.2 

Flow  Curve  for  the  "Rigid"  Pseudo-Plastic  Material 


back-pull,  and  hydraulic  extrusion  with  zero  drawing  force.  Numerical 
results  are  presented  for  hydraulic  extrusion  and  for  wire  drawing,  with 
n  =  0.1;  only  the  case  of  no  back-pull  is  presented,  since  computations 
for  the  case  including  back-pull  are  excessively  laborious. 

The  wire  drawing  analysis  shall  now  be  presented.  FIGURE  4.3 
serves  to  clarify  the  physical  situation  and  to  introduce  the  necessary 
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notation;  it  is  further  noted  that 


n  i  2/3* 


in  the  following  analysis. 


*The  case  n  =  2/3  is  excluded  since  the  relation  (4.20),  which  is  used 
extensively  in  the  numerical  procedure,  does  not  hold  if  n  =  2/3. 


u  ar  rtorriw  t(0S.£)  norJelsn  son>a  bsbufDxs  sr  £\S  3  n  9263  sdT* 


68 


4.1  Wire  Drawing 

The  entry  and  exit  sections,  AB  and  CD  respectively,  as  shown 
in  FIG.  4.3,  are  taken  as  plane  sections,  rather  than  circular  arcs;  this 
requires  a  knowledge  of  a  .  The  tensor  transformation  rule  from  spherical 
polars  to  cartesian  coordinates  gives 

aZ=(~Hr^)  +  cos  2(^  ‘  Trcj)Sin  2<*  *  (4-2) 


The  polar  components  of  the  stress  tensor  are  calculated  from 

a.  .  ~  s .  .  —  p6 •  •  . 

1J  1J  1J 

The  deviatoric  and  isotropic  components  are  given  by  equations  (3.6)  and 
equation  (3.8)  respectively,  for  completeness  these  equations  are  reproduced 
below 

5-n 

0  2  An 

sr  =  3nU  mG  ’ 
r 

3-n 

s4>  =  se  =  r3npA  mG  *  ^3-6^ 

1-n 

.  .  2  2  uAn  ,r 

sr<t>  "  r3n  G  ’ 


and 


1-n 

2  p 

p  =  2-3nliA-  [-(l-n)mG  -  {(m'G)1  +  m'Gcot^}  ]  +  C  ,  (3.8) 


3n 
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where 


n-1 


G  =  [12m2  +  (m1)2] 


In  expanded  form  the  stress  tensor  components  become 


1  -n 

-  2  2  yAn  [— {-  —  mG  +  ^-  ( (m'G) '  +  m'Gcotcf))}  -  H]  , 

1-n 

■  2  2  uAn  2(3„~2)  mG  +  ^rUm'G)'  +  m'Gcot*)}  -  H]  , 


=  CTe 


3n 


(4.3) 


1  -n 

Tr<(>  ■  2  2  <*"  (m'G)] 


The  constant  C  in  equation  (3.8)  has  been  replaced  by 


1-n 

H  2  2  yAn  ; 


the  constant  H  is  determined  by  the  back-pull  condition. 

Substitution  of  equations  (4.3)  in  equation  (4.2)  yields 


o_  =  2 


-  H  ]  . 


(4.4) 


where 


f(<t)  =  3mG(  1  -cos2«J>)  -  ^  {(m'G) '  +  m'Gcot<t>}  -  m'Gsin2<J>  .  (4.5) 
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If  the  mean  tensile  stress  over  the  entry  section  AB  is  given  by  Tq  , 
where  T0  must  be  less  than  the  tensile  yield  stress,  then  the  condition 
of  back-pull  becomes 


whence 


d(area) 

at  AB 


» 


H  = 


z03ntan2a 


•a 


f (<f>)  cos3n“2cf>  tancj)d<t>  - 


1-n 

,  2 


(4.6) 


UA 


n 


In  equation  (4.6),  zQ  is  the  distance  from  the  virtual  apex  of  the  cone, 

0,  to  the  entry  section  AB,  as  shown  in  FIG.  4.3.  Finally,  from  equations 
(4.4)  and  (4.6),  the  expression  for  a in  terms  of  TQ  is 


“,-2 


z0^ntan^a 


■a 


f(<f>)  cos3n‘2cl)tan(J)d(j)]  +  Tq  .  (4.7) 


The  basic  unknown  of  the  problem  is  the  drawing  stress  T-|  ,  which 
must  be  less  than  the  tensile  yield  stress.  It  is  found  from  the  relation 


T1  = 


ttD 


1  CD 


at  CD 


d(area)  , 


which  yields 


3-n 
,  2 


Ti  ■ 


uA 
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z-j  3ntan2a 


z,  3n  ra 
[1  -  (y1)  ] 
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f (4> )  cos3n”2(})tancj)d(|)  +  T 
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where  z-j  is  the  distance  from  the  virtual  apex  of  the  cone  to  the  exit 
section  CD,  as  shown  in  FIG.  4.3.  In  terms  of  the  diameters  DQ  and  D-j 
of  the  entry  and  exit  sections  respectively,  the  expression  for  T-j  becomes 


3+5n 


T-,  = 


nn.  3n-2 
yA  tan  a 


D 


3n 


[i  - 


D-,  3n 

$  i 


ra 


f(<f>)  cos ^d>tan4>d4>  +  T 


Upon  the  introduction  of  the  reduction  of  the  die,  R,  given  by 


initial  area  -  final  area 

initial  area 


■  [i  -  (q1)  ]  • 
o 


gives 


3+5n 


T,  = 


nn,  3n-2 
yA  tan  a 

pn  3n 


3n 

[1  -  (l-R)2  ] 


■a 


f(<i>)  cos^n"^4)tan(j)d4)  + 


(4.8) 


The  drawing  stress,  T-| ,  may  not  exceed  the  tensile  yield  stress,  and  its 
maximum  value  will  now  be  calculated. 

As  previously  mentioned,  the  minimum  value  of  I  occurs  at  the 
point  E,  see  FIG.  4.3,  and  is 


I  = 


2/3  A 


in  terms  of  the  diameter  D-j  ,  and  the  reduction  R,  this  becomes 
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16/3  Atan  a 

D3 


(1-R) 


3/2 


(4.9) 


To  obtain  the  critical  value  JQ,  the  equation  linking  the  in¬ 
variants  J  and  I  is  required;  this  is  obtained  by  squaring  both  sides  of 
equation  (1.6),  and  yields 


1-n 

J  *  2  2  pln  . 


(4.10) 


Substitution  of  equation  (4.9)  in  equation  (4.10)  results  in 

7n+l  n 

n  "?n 

3n/2 


0o  =  1  2  -..323nMAntan3na  (1.R) 


D 


1 


As  previously  noted  in  equation  (4.1),  the  maximum  tensile  yield  stress  is 


Y  =  /3  J 

o 


» 


which  becomes 


7n+l  n+1 


Y  = 


0  2  n  2  An.  3n 
2  3  uA  tan  a 


D 


3n 


(1-R) 


3n 

2 
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Since  T-|  may  not  exceed  Y,  the  following  inequality 
3+5n 


3n 


9  2  AnH-=«3n-2  o~ 

2 - yAJan - .  (1_R)2 


D 


1 


■a 


f (<J>)  cos^n"^(j)tan(i)d4)  +  T 


o 


o 


7n+l  n+1 


o  2  0  2  An. 3n  0 

<  2  3  yA  tan  a  ^  _  ^2 
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3n 


(4.11) 
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determines  the  maximum  reduction  possible;  apart  from  any  negative  normal 
stresses  at  the  walls,  which  may  occur  during  the  process .  For  zero 
back-pull,  from  equation  (4.11),  for  the  maximum  reduction: 

n+1  _2 

R  <  1  -  [ii11-1-3-— ).tJ.022 -  +  1  ]  3n  . 

/a  f(<|>)  cos^n"^0tancj)dct) 
o 

An  expression  for  the  die  pressure,  P^.  ,  shall  now  be  developed; 
this  expression  could  be  useful  in  a  stress  analysis  of  the  die  proper. 

By  definition 


P 


die 


0 


The  stress  component  at  an  arbitrary  point  P  is  obtained  by  substitution 
of  equation  (4.6)  in  the  second  of  the  set  of  equations  (4.3),  which  re¬ 
sults  in 


1-n 

0  2  An  r  1 
Qcj)  =  a0  =  2  pA 

2 

"  3n ,  2 

z  tan  a 
o 

Upon  substitution  of  <j>  =  a 


{2_(3n-?]_  mG  +  1_  ( (m'G) '  +  m'Gcot<j>)} 
^  'In-? 

f (<f>)  cos  q>tan<t>d4)]  +  TQ  . 
o 

in  equation  (4.12),  and  noting  that 


(4.12) 


m(a)  =  0  , 
n-1  j 

4>=a 


G(a)  =  (m‘) 


9 


. 


74 


G 1 (a)  =  (n-1 ) (m ' )n  2m 


i  i 


<|>=a 


the  die  pressure  becomes 


1-n 

Pdie  =  2  2  vA^'~k  {3F  (n(m')n'1  m"  +  (m'  )ncot<j>)  } 

s  <j)sa 


(4.13) 


z  3ntan2a 
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f ( 4>)  cos3""2  4>tan4)d(j)]  -  T 


The  parameter  s  is  measured  along  OCA,  as  shown  in  Fig.  4.3,  and  is  calcu¬ 
lated  from 


s  =  Y 


1 


cosa  * 


(4.14) 


where  y  is  a  variable  with  the  range  of  values  given  by 


z  -1/2 

1  <  y  <  —  =  (1-R) 


—  '  —  z 


1 


Under  the  substitution  (4.14)  the  expression  for  P^.  ,  in  terms  of  the 
reduction  R,  becomes 


l+5n 
2  nn 


n 


die 


^  <4-  (n(m')'1'1  m"  +  (m' )"cot<j>)  } 
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3n  3n 
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+  2(1-R)2  tan3n'2  a 


a 


o 


4)=ot 


f(<J>)  cos3""2  <j>tan<j>d<j>]  -  TQ  . 


(4.15) 


Similarly,  the  expression  for  the  shear  stress  at  the  boundary  is  given  by 
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l+5n 

TrJ  °  3nyA3n  [(m'  )nsin3nj)]  .  (4.16) 

$  <(>=a  yJn  <(>=cx 

It  is  noted  that  the  integral 

ra 

f (<t>)  cos  4>tan4>d(j) 
o 


occurs  in  many  of  the  preceeding  expressions.  It  is  not  necessary  to  per 
form  this  integration  numerically  because  the  requirement  of  static  equili 
brium  of  the  material  between  the  entry  and  exit  sections  yields  an  ex¬ 
pression  for  this  integral  in  terms  of  a,  m'(a)  ,  and  m''(a)  .  This  ex¬ 
pression  is  now  obtained.  The  net  axial  component  of  the  die  force,  F^.  , 
is 


F 


die 


z  /cosa 

9  f  0 

2-rrsin  a  j  sds  , 

s=z-j  /cosa 


(4.17) 


which,  upon  substitution  of  equation  ( 4 „ 1 3 ) ,  reduces  to 
5n-l  3n-2 
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The  net  axial  component  of  the  shear  force,  F 
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F 


shear 


zQ/cosa 

2tts  i  nacosa 


<|>=a 


sds 


9 


s-z^/cosa 


which,  upon  substitution  of  t  , |  » reduces  to 

^  cf>=a 

5n-1  3n-2  ] 

Fshear  =  2  2  ^  Dl2'3TO-R)  3  -  l][(m')n  -■/■■2.3nfCOS<t']  -(4.18) 

v  '  cj)=a 

The  net  drawing  force,  Fdrawing  ,  is 

itD12T1 

^drawing  4  * 


which,  upon  substitution  of  the  expression  for  T-j  ,  equation  (4.8),  reduces 
to 


drawing 
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5n-1  3n 
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The  condition  of  static  equilibrium  becomes 


ttD  T 

F  +  F  +  — - — -  F  =0 

shear  die  4  drawing 


(4.19) 


Upon  substitution  of  the  above  obtained  expressions,  there  results  from 
equation  (4.19) 
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(4.20) 


It  is  sometimes  useful  to  treat  the  boundary  effects  by  the 
introduction  of  a  coefficient  of  friction  as  is  done  for  Coulomb  friction; 
this  has  been  done  in  various  manners  by  Sachs  (34),  and  also  by  Shield  (31). 
Proceeding  in  the  fashion  of  Shield,  an  average  coefficient  of  friction, 
defined  by 
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which  becomes  in  expanded  form 
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Results  for  wire  drawing  are  presented  in  FIGS.  6.A-1  and  6.A-2. 
As  previously  mentioned,  TQ  was  set  equal  to  zero  to  facilitate  the 
numerical  computations. 

4.2  Hydraulic  Extrusion 

The  hydraulic  extrusion  of  the  pseudo-plastic  material,  under 
conditions  of  zero  drawing  force, that  is 


CD 


at  CD 


d(area)  =  0  , 


(4.22) 


is  now  considered.  From  the  above  expression  (4.22)  the  constant  H  in 
equation  (4.6)  is  found  to  be 


H  =  — f  f(4>)  cos3n"2<l>tan<J>d4). 


z-,3ntan2a  * 
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(4.23) 


In  terms  of  H, 
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where  f(cj>)  is  defined  by  equation  (4.5). 

The  basic  unknown  for  the  hydraulic  extrusion  is  the  ram  pressure 
Pram*  this  compares  to  T-|  of  the  drawing  process. 

The  equilibrium  condition  at  the  entry  section  AB  is 
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from  which  there  is  obtained  in  terms  of  the  reduction  R 
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Consequently,  the  ram  force,  defined  by 
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From  the  second  of  the  set  of  equations  (4.3)  and  using  equation 
(4.23),  and  evaluating  at  <f>  =  a,  there  results 
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The  net  axial  component  of  the  die  force,  obtained  by  substituting  equation 
(4.24)  in  equation  (4.17)  is 
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The  expression  for  t 
(4.16) ,  and  is 

l+5n 


rcj>  a 
cj)=a 


is  as  previously  obtained  in  equation 


rcj) 


2  dA  r/«.i\ir'-.*w.3nin 

3n  5n  '  sin  ^ 


4>=a  D-j  y 


(j)=a 


(4.16) 


The  net  axial  component  of  the  shear  force,  therefore,  remains  unchanged 
and  is  given  by  equation  (4.18),  which  for  completeness  is  reproduced; 
thus , 

5n-l  3n-2 

Fshear  =  2  2  ^  Dl2"3"  [{1'R)  2  ‘  •  <4-18) 

It  is  noted  that  the  overall  statical  equilibrium  condition  given  by 

equation  (4.19)  is  essentially  independent  of  TQ  and  T^ ,  and  therefore 

of  P  ,  ,  hence  the  integral 
ram  J 


a 


f((j))  cos^n"^cj)tan(j)d(}) 


is  defined  in  terms  of  boundary  values  by  equation  (4.20). 

The  average  coefficient  of  friction,  faverage»  as  defined  by 
equation  (4.21)  becomes  in  expanded  form 
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FIG.  6.A-3shows  the  relationship  between  P  ,  and  semi-angle  a  for  various 

r  ram  3 

reductions . 

It  is  important  to  note  that  all  calculations  for  both  wire 
drawing  and  hydraulic  extrusion  are  based  on  a  constant  diameter  at  the 
exit  section. 


CHAPTER  V 


EXPERIMENTAL  INVESTIGATION 

The  experimental  investigation  of  the  axially  symmetric  con¬ 
verging  flow  of  certain  non-Newtonian  materials  in  a  circular  cone  is 
presented  in  this  chapter. 

The  purpose  of  this  experimental  work  is  to  investigate  the 
kinematic  behavior  of  creeping  flow  of  non-Newtonian  materials.  This 
investigation  is  concerned  mainly  with 

(i)  determining  if  the  creeping  flow  of  non-Newtonian  materails  in 
a  circular  cone  yields  a  radial  flow  pattern,  and 
(ii)  observing  whether  or  not  the  material  adheres  to  the  walls  of 
the  die  during  the  extrusion  process. 

The  assumption  that  the  relative  motion  of  a  non-Newtonian  material 
in  contact  with  a  solid  body  is  zero,  may  not  always  be  justified;  this 
becomes  more  evident  if  pure  lead  is  considered. 

Often  pure  lead  is  used  in  experiments  whose  purpose  is  to 
examine  the  validity  of  solutions  for  a  rigid  plastic  solid.  At  a  perfectly 
rough  boundary  the  solutions  for  flow  of  a  rigid  perfectly  plastic  solid 
indicate  that  there  can  be  slip  at  the  boundary,  that  is  the  relative 
motion  is  non-zero.  The  problem  of  flow  in  a  converging  channel  (30)  is 
an  example  of  this.  In  this  solution,  as  in  other  rigid  plastic  solutions 
where  there  is  relative  motion  between  a  perfectly  rough  boundary  and  the 
material,  the  shear  strain  rate  is  infinite.  This  is  permissible  in  an 
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inviscid  material  but  is  not  permissible  for  a  material  such  as  a  Bingham 
material  that  exhibits  viscous  effects.  The  classical  perfectly  plastic 
solid  is  not  rate  dependent.  It  has  been  shown  by  Luming  (35)  that  at 
room  temperature  the  behavior  of  pure  lead  can  closely  be  approximated 
by  that  of  a  Bingham  solid,  with  a  low  coefficient  of  viscosity.  Conse¬ 
quently  pure  lead  may  be  regarded  as  a  non-Newtonian  material.  Possibly, 
the  amount  of  slippage  at  the  boundary  should  be  considered  as  a  function 
of  the  Bingham  number*;  this  function  would  be  of  such  a  nature  that  the 
amount  of  slippage  at  the  walls  increases  with  increasing  Bingham  number. 

In  this  thesis  these  considerations  have  not  been  taken  into  account,  and 
the  classical  assumption  of  no-slip  at  a  boundary  is  used.  Further  experi¬ 
mental  work  is  needed  to  study  the  behavior  of  non-Newtonian  materials 
in  contact  with  solid  boundaries. 

5.1  Experimental  Apparatus 

The  experimental  investigation  consisted  essentially  of  extruding 
lead  and  wax  specimens.  A  rigid  conical  die,  of  circular  cross-section, 
was  designed  to  perform  these  experiments.  The  extruding  apparatus  con¬ 
sisted  of  three  main  parts: 

(i)  The  Rigid  Die  and  Separation  Plate-  the  die  was  designed  so  that 

(a)  removal  of  the  specimens  was  possible,  and 

(b)  it  could  be  used  to  prepare  the  specimens,  so  that  the 
specimens  were  in  two  halves  with  a  grid  showing  on  the 

*The  Bingham  number  is  defined  as  B  =  kL/yll,  where  k,L,y,  and  U  are  the 
yield  stress  in  pure  shear,  the  characteristic  length,  the  Newtonian 
coefficient  of  viscosity,  and  the  characteristic  velocity,  respectively  (36). 
The  Bingham  number  is  a  measure  of  the  ratio  of  the  relative  magnitudes 
of  the  plastic  forces  to  the  viscous  forces. 
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face  of  one  of  the  halves,  thus  facilitating  in¬ 
vestigation  of  the  flow  behavior. 

In  order  to  facilitate  the  removal  of  the  specimens,  the  die,  approxi¬ 
mately  8"  long,  was  of  the  split  type,  bolted  together  with  ten  7/16" 

N.C.,  and  two  3/8"  N.C.  bolts,  as  shown  in  FIG.  5.1.  To  ensure  perfect 
alignment  of  the  two  halves,  the  halves  were  provided  with  some  1/4" 
aligning  pins.  The  semi-angle  of  the  circular  cone  was  5°,  with  an  exit 
diameter  of  1/4".  Due  to  the  small  exit  diameter,  the  machining  of  the 
cone  had  to  be  done  in  two  stages;  therefore,  the  die  had  to  be  built  up 
out  of  two  parts  as  is  shown  in  FIG.  5.1. 

To  prepare  the  specimens,  a  1/4"  plate  was  prepared  to  be  clamped 
between  the  halves  of  the  die;  this  way  two  molds  were  formed,  so  that  the 
specimen  halves  could  be  poured  with  liquid  lead  or  wax.  Furthermore,  the 
plate,  provided  with  the  appropriate  holes  for  the  aligning  pins,  was  pre¬ 
pared  with  a  grid  on  one  of  its  sides;  the  grid  consisted  of  radial  lines 
at  intervals  of  1°  -  40',  and  with  circular  arcs  about  the  virtual  apex 
of  the  cone  at  1"  intervals,  see  FIG.  5.2.  A  circular  ring,  of  semi¬ 
circular  cross-section,  was  made  in  the  top  of  the  die;  this  provided  the 
location  for  the  oil  ring,  which  was  clamped  between  the  die  and  the  piston 
housing ,  see  FIG.  5.1. 

(ii)  The  Piston  Housing  -  the  piston  housing  consisted  essentially 
of  a  thick  cylinder,  approximately  3"  long,  see  FIG.  5.3.  The  end  of  the 
piston  housing,  which  is  to  be  bolted  to  the  die  with  six  7/16"  N.C.  bolts, 
was  provided  with  a  groove  similar  to  the  one  on  the  top  of  the  die.  The 
outside  diameter  was  the  same  as  that  of  the  die,  that  is  5".  The  inside 
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6-1/2" 


1-1/2" 


FIGURE  5.1 

Typical  Half  of  the  Die 
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diameter,  for  the  piston  to  be  fitted  in,  was  machined  at  1  7/16". 

(iii)  The  Piston  -  the  piston  consisted  mainly  of  a  piece  of  shaft, 

5"  long,  of  circular  cross-section  with  a  diameter  of  1  3/8".  At  one  end 
of  the  piston  four  grooves  were  provided  to  house  the  oil  rings,  as  shown 
in  FIG.  5.3.  A  hole  was  bored  in  the  piston,  and  a  hydraulic  pump  was  at¬ 
tached;  this  enabled  the  extrusions  to  take  place  continuously  without  the 


Piston 


To  Hydraulic 
Oil  Pump 


Passage  for  Oi 1 
Addi tions 


Oil  Rings 


Holes  to  Take 
7/16"  N.C.  for 
Mounting  on  Die 


Piston  Housing 


FIGURE  5.3 


The  Piston  and  the  Piston  Housing 
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disadvantage  of  removing  the  piston  in  order  to  fill  up  the  oil  level 
during  the  extrusion. 

5.2  Results  of  Experiments 

As  previously  mentioned,  the  materials  used  in  this  investi¬ 
gation  were  lead  and  wax.  The  lead  used  was  that  as  supplied  by  Canada 
Metals  Ltd.  Its  assay  is  0.02%  antimony,  0.015%  arsenic,  and  99.965%  lead. 
The  material  constants  were  found  experimentally  by  Luming  (35)  as 
p  =  1372  psi-sec.,  and  k  =  1170  psi . ,  where  y  is  the  Newtonian  coefficient 
of  viscosity,  and  k  is  the  yield  stress  in  pure  shear.  The  waxes  used 
were 

(i)  Parowax  -  made  by  the  Imperial  Esso  Oil  Company  -  of  which 
the  material  constants  were  not  known,  and 
(ii)  Bee's  wax,  similarly  without  a  knowledge  of  its  properties. 

Before  discussing  the  actual  extrusions,  the  preparation  of  the  specimens 
is  considered  briefly. 

The  Specimens 

The  material,  lead  or  wax,  was  melted  and  poured  into  the  molds; 
as  previously  discussed,  these  molds  were  obtained  by  the  insertion  of  the 
separation  plate  between  the  two  die  halves.  The  formation  of  pipes, 
longitudinal  cavities  due  to  shrinkage  toward  the  side  of  the  mold  during 
solidification,  was  minimized  by  adding  molten  material  into  the  cavities 
during  the  cooling  process.  This  technique  was  used  throughout,  both  for 
the  wax  and  the  lead.  After  removing  the  die  halves  from  the  separation 
plate,  some  specimens  had  to  be  abandoned,  in  particular  the  waxes,  since 
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sticking  of  the  material  to  the  separation  plate  caused  some  specimens  to 
crumble. 

Lead  Extrusions 

After  the  die  was  assembled,  the  two  halves  of  the  specimen  were 
tapped  in  simultaneously,  so  that  a  proper  fit  could  be  ensured.  The 
apparatus  was  completed  by  the  mounting  of  the  piston  housing,  with  the 
appropriate  oil  ring  inserted.  With  the  insertion  of  the  piston,  extru¬ 
sion  could  now  take  place.  The  applied  load  increased  to  approximately 
4500  lbs.,  at  which  point  the  oil  ring,  separating  the  die  and  the  piston 
housing,  could  not  support  a  further  increase  in  load,  consequently  the 
oil  was  seeping  from  the  sides  and  extrusion  of  the  specimen  did  not  take 
place.  A  few  more  attempts  were  made,  but  all  with  the  same  result.  Before 
attempting  any  more  extrusions, the  piston  housing  was  milled  down  0.003" 
so  that  the  piston  housing  could  be  bolted  to  the  die  at  a  higher  torque, 
thus  providing  a  tighter  fit  for  the  oil  ring.  Extrusions  were  attempted 
again,  but  the  problem  of  oil  seepage  remained.  Before  abandoning  the  lead 
extrusions  completely,  a  lead  plug  was  machined  to  fit  between  the  specimen 
and  the  piston;  this  plug  was  used  instead  of  the  hydraulic  fluid,  thus  re¬ 
moving  the  sealing  problems  of  the  previous  extrusions.  After  the  apparatus 
was  assembled,  extrusion  could  take  place  again.  The  applied  load  increased 
to  approximately  45,000  lbs.  before  plastic  flow  occurred.  Although  some 
minor  flashing  occurred,  due  to  a  small  separation  of  the  die  halves,  the 
specimen  was  kept  for  investigation  of  the  flow  behavior. 
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Wax  Extrusions 

The  hydraulic  extrusions  using  Parowax  were  attempted  first. 
Before  extrusion,  some  coloring,  using  a  felt  pen,  was  applied  to  the  grid 
on  the  specimen  so  that  after  extrusion  the  grid  would  be  more  distinct, 
thus  providing  easier  evaluation  of  the  flow  behavior.  The  applied  load 
increased  to  approximately  4000  -  5000  lbs.,  and  plastic  flow  did  indeed 
occur,  however,  the  extrusion  was  not  as  expected.  Due  to  the  contraction 
of  the  specimen  halves  on  cooling,  the  flat  sides  in  contact  with  each 
other  during  extrusion  did  not  seal  properly,  consequently  the  oil  was 
forced  through  this  opening,  and  extruded  the  bottom  portion  of  one  of  the 
halves  of  the  specimen.  Further  attempts  were  made,  but  all  with  the 
same  result. 

In  order  to  overcome  this  difficulty,  the  specimens  were  poured 
to  approximately  1/2"  from  the  top,  and  after  the  two  halves  were  fitted 
in  the  die,  a  wax  plug  was  poured  on  top  of  the  specimen,  thus  sealing  any 
imperfections  which  caused  oil  seepage  problems  during  the  previous  ex¬ 
trusions.  Extrusions  took  place  again  with  only  a  few  useful  results.  It 
appeared  that  somehow  the  oil  bypassed  the  plug  causing  either  the  failures 
as  described  above,  or  due  to  the  uneven  pressure  on  the  specimen  halves, 
parts  of  the  specimen  "cork-screwed"  out  leaving  the  other  half  crumbled 
in  the  die. 

Most  of  the  few  extrusions  which  did  indeed  give  good  plastic 
flow  were  of  no  value,  since  the  two  specimen  halves  were  fused  together, 
and  due  to  its  brittle  nature  could  not  be  separated  wholly.  A  few  speci¬ 
mens  were  kept,  however,  and  investigated. 
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The  extrusions  using  bee's  wax  gave  the  same  difficulties  as 
the  ones  with  Parowax,  except  not  to  such  a  great  extent.  Numerous  ex¬ 
trusions  were  attempted,  all  without  success.  Another  complication,  which 
was  not  the  case  with  the  Parowax,  was  that  the  coloring  of  the  specimen 
with  a  felt  pen  was  not  effective.  The  coloring  of  the  specimen  was 
changed  to  the  application  of  graphite  in  the  grid  on  the  separation  plate 
before  pouring  the  specimens;  this  indeed  resulted  in  a  good  grid.  More 
attempts  were  made  to  obtain  good  plastic  flow  of  the  whole  specimen,  most 
without  good  results;  however  one  of  such  extrusions  did  result  in  a  speci¬ 
men  which  was  used  for  the  investigation  of  the  flow  behavior. 

5.3  Conclusions 

The  lead  specimen,  which  was  retained  for  investigation,  was 
obtained  from  the  extrusion  during  which  some  minor  flashing  occurred.  The 
specimen  was  investigated,  and  it  was  observed  that  the  radial  lines  re¬ 
mained  radial.  Also  the  presence  of  relative  motion  between  the  specimen 
and  the  die  was  observed,  notwithstanding  the  fact  that  the  walls  of  the 
die  were  degreased  before  extrusion. 

The  specimens  obtained  using  Parowax  were  investigated,  and  it 
was  found  that  radial  lines  did  indeed  remain  radial,  thus  predicting  a 
radial  flow  pattern,  at  least  for  creeping  flow,,  Furthermore,  it  was  ob¬ 
served  that  the  material  did  adhere  to  the  walls  of  the  die,  as  is  the 
classical  assumption  for  Newtonian  viscous  flow.  The  same  flow  behavior 
was  observed  with  the  specimen  obtained  from  the  bee's  wax  extrusion. 

Since  the  specimens  of  the  Parowax  were  obtained  during  the 
earlier  stages  of  the  experimental  period,  no  photographs  were  taken  since 
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it  was  assumed  that  more  and  better  specimens  would  be  obtained  once  the 
technique  was  ironed  out.  Furthermore,  used  specimens  were  melted  and 
used  over  again.  Since  no  other  successful  extrusions  were  attained 
during  the  many  experiments  performed  thereafter,  no  photographs  are 
available  to  include  in  this  chapter.  The  extrusions  using  bee's  wax, 
which  were  performed  during  the  final  stages  of  the  experimental  period, 
did  yield  one  specimen  which  extruded  rather  well.  However,  the  two 
halves  fused  together,  therefore,  one  of  the  halves  was  cut  off,  taking 
care  not  to  damage  the  grid.  The  grid  was  clear  enough  for  investigation, 
and  showed  both  radial  flow  and  adhesion  to  the  walls  of  the  die.  Al¬ 
though  the  visibility  of  the  grid  was  good,  the  photographing  of  the 
same  was  not  clear  at  all,  and  the  photographs  obtained  do  not  merit  in¬ 
clusion  in  this  chapter. 

It  may  be  concluded  from  this  experimental  investigation,  there¬ 
fore,  that  the  flow  pattern  for  the  creeping  flow  of  the  non-Newtonian 
materials  tested  is  radial,  at  least  for  small  semi-angles.  This  con¬ 
clusion  is  not  necessarily  valid  for  larger  semi -angles;  more  extensive 
experimental  work  is  required  to  investigate  any  possible  deviations  from 
radial  flow,  which  may  result  from  larger  semi-angles. 

Further  experimental  investigations  of  the  creeping  flow  of 
lead,  which  may  be  assumed  to  behave  as  a  Bingham  solid  for  a  certain 
range  of  strain  rates,  would  be  a  significant  contribution,  especially 
if  the  theoretical  investigation  is  considered  in  greater  detail. 
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CHAPTER  VI 


CONCLUDING  REMARKS 

The  main  purpose  of  this  chapter  is  to  consider  further  the 
pseudo-plastic  flow  problems,  which  were  discussed  previously  in  this 
thesis. 

The  numerical  results  of  these  flow  problems  are  incorporated 
in  the  first  section  of  this  chapter,  while  the  results  of  wire  drawing 
and  hydraulic  extrusion  are  presented  in  the  section  dealing  with  those 
processes.  A  discussion  on  how  closely  these  results  approximate  the 
correct  solution  is  given.  Furthermore,  the  usefulness  of  these  results 
is  investigated  by  consideration  of  the  objections,  raised  by  Reiner  (24), 
to  the  use  of  the  power  law  model,  which  is  used  in  this  thesis.  The 
results  of  the  application  of  the  axially  symmetric  pseudo-plastic  flow 
to  wire  drawing  and  hydraulic  extrusion  are  then  briefly  considered.  In 
conclusion,  a  method  due  to  Sutterby  (37),  which  could  possibly  give  an 
approximate  solution  to  the  axially  symmetric  flow  of  a  Bingham  solid 
in  a  circular  cone  of  small  semi -angle,  and  of  finite  length,  is  briefly 
discussed. 

6.1  Results 

Various  results  of  the  numerical  procedure  have  been  presented 
at  the  end  of  this  section.  The  results  for  the  velocity  profiles,  with 
semi-angles  15°,  30°,  45°,  and  60°,  and  for  various  values  of  n,  are  pre¬ 
sented  in  FIGS.  6.1 1-2  to  6,11-5  and  in  FIGS.  6.III-2  to  6.III-5  for  the 
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plane  and  axially  symmetric  flow,  respectively.  It  is  recalled  that  a 
knowledge  of  m ' 1 ( 0 )  was  required  to  solve  the  boundary  value  problems. 

The  relationships  between  this  parameter  and  n  for  various  semi -angles  a 
are  presented,  therefore,  in  FIG.  6.II-1  for  the  plane  flow,  and  for  the 
axially  symmetric  flow  in  FIG.  6.III-1. 

The  velocity  profiles  for  the  limiting  cases,  n  =  0,  and  n  =  1, 
are  shown  for  comparison.  The  solution  for  the  case  n  -  1,  corresponding 
to  the  Newtonian  viscous  liquid,  may  be  written  as 

2  2 

/  _  cos  <t>  -  cos  a 

m  \  (p )  2  • 

sin  a 

which  is  according  to  the  expression  originally  obtained  by  Harrison(38) 
and  Bond  (39).  For  the  plane  flow,  it  can  be  shown  that  the  solution 

2  2 

_/A\  cos  0  -  cos  a 

m(0)  =  - 2 - 

sin  a 

satisfies  the  governing  differential  equation.  The  solution  for  the  case 
n  =  0,  corresponding  to  the  perfectly  plastic  von  Mises  solid,  for  the 
plane  flow  case  may  be  written  as 

”  F^coTTijr 

where  ip  is  the  angle  between  a  radius  and  the  direction  of  the  algebraical ly 
greatest  principal  stress,  and  the  parameter  c  is  defined  by 
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..  c  ......  tan*1  /P-l  =  i  it  +  a 

,£z-7  FTl  4 

for  the  case  of  perfectly  rough  walls;  these  equations  are  due  to  Hill  (30). 
For  the  axially  symmetric  case,  with  n  =  0,  the  velocity  field  was  ob¬ 
tained  only  partially,  since  the  derivatives  at  the  walls  are  infinite. 

The  value  used  for  m ' ' ( 0 ) ,  needed  to  start  the  integration  procedure, 
was  that  obtained  by  extrapolating  the  curves  in  FIG.  6.III-1  to  n  =  0. 

It  is  noted  from  the  velocity  profiles  that  there  is  a  distinct 
difference  between  the  curves  n  =  0.1  and  n  =  0.  This  difference  is  more 
pronounced  than  that  between  the  curves  n  =  0.9  and  n  =  1.  To  illustrate 
this  further,  the  velocity  profiles  for  the  plane  flow,  with  n  =  0.04 
and  n  =  0.07  have  been  presented  for  a  channel  with  a  semi-angle  of  15°, 
see  FIG.  6. 1 1-2.  This  marked  difference  can  be  attributed  to  the  fact 
that  in  the  range 


0  <  n  <  1 


the  material  is  assumed  to  adhere  to  the  walls,  whereas  for  the  case 
n  =  0,  there  is  slippage. 

In  the  solutions  of  the  stress  fields  for  the  plane  perfectly 
plastic  von  Mises  flow  in  a  converging  channel  due  to  Nadai  (29),  and  that 
for  the  axially  symmetric  flow  in  a  circular  cone,  due  to  Shield  (31),  it 
was  assumed  that  for  perfectly  rough  walls,  the  shear  stress  at  the  walls 
is  k,  where  k  is  the  yield  stress  in  pure  shear.  This,  however,  is  not 
the  case  if  a  viscous  medium  is  considered;  the  shear  stress  at  the  wall 
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is  determined  by  the  velocity  field.  To  illustrate  this  further,  the 
results  presented  in  FIGS.  6,11-6  and  6.III-6  for  the  plane  and  axially 
symmetric  flow,  respectively,  show  that  for  n  =  0  the  shear  stress  at 
the  wall  is  independent  of  the  semi -angle  a,  whereas  for  values  of  n 
in  the  range 


0  <  n  <  1  , 


the  shear  stress  at  the  wall  is  dependent  on  the  semi-angle  of  the  cone. 
Furthermore,  this  dependence  increases  with  increasing  n.  For  complete¬ 
ness,  results  of  the  relationships  between  the  shear  stress  at  the  wall 
and  semi-angle  a  have  been  presented  in  FIGS.  6.II-7  to  6,11-9  for  the 
plane  flow,  and  in  FIGS.  6. 1 1 1-7  to  6. 1 1 1-9  for  the  axially  symmetric 
flow.  It  is  noted  that  in  FIGS.  6.II-9  and  6,111-9  the  non-dimensional 

2n  3^ 

quantities  xrQ(a)  r  /yQr  and  Tr^(a)  r  /yQn  have  been  replaced  by 
2  3  n  n 

ax^a)  r  /yQn  and  ax^fa)  r  /yQn  for  the  plane  and  axially  symmetric 
flow,  respectively.  This  replacement  is  necessary,  since  the  results 
of  x^ta)  r  /yQn  versus  a  for  the  plane  flow,  and  those  of  i^fa)  r"n/ 
yQn  versus  a  for  the  axial ly -symmetric  flow  yield  curves  which  are  of 
such  a  nature  that  only  a  qualitative  and  not  a  quantitative  study  can 
be  made. 


Since  the  solutions  to  the  axially  symmetric  flow  problems, 
with  n  =  0.1,  are  applied  to  wire  drawing  and  hydraulic  extrusion,  the 
stress  distributions  in  the  material  are  presented  only  for  that  parti¬ 
cular  value  of  n.  The  semi-angles  for  this  purpose  were  chosen  to  be 


' 


97 


15°,  30°,  45°,  and  60°.  Although  large  semi-angles  are  of  no  practical 
consequence  in  the  application  to  wire  drawing,  they  are,  nevertheless, 
of  academic  interest,  and  thus  merit  presentation.  These  results  are 
presented  in  FIGS.  6 . 1 1 1 - 1 0  to  6.III-13.  Similar  results  were  obtained 
for  the  plane  flow  case,  and  are  presented  in  FIGS.  6.II-10  to  6. 11-13. 


tst.19.-  (  tvAi  bn.  .Jesnstnf  0'msbt.m  10 

,,rtr^  afiftTn  fO< 


■ 


RESULTS 


OF 

PLANE  PSEUDO-PLASTIC  FLOW 


IN  A  CONVERGING  CHANNEL 


30 

25 

20 

15 

10 

5 

0 


99 


_ 

c 

1 

oC  =  6C 

0 

o 

0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 

n 

FIGURE  6.11-1 


m"(o)  versus  n 


m  (e) 


100 


e 

FIGURE  6.E-2 
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FIGURE  6.H-3 
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FIGURE  6.H-5 
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FIGURE  6.31-6 

Tro  [oC )  r2n/  ju  Qn  versus  n 
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FIGURE  6.1-8 

Tr0  (06)  r 2n/ ju  Qn  versus  oi 
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FIGURE  6.IC-9 

o£Tre  (o£)  r2n/  juQn  versus 


o L 


/  juQ 


4 


108 


FIGURE  6  .  It-  10 

CXjj  r2n/juQn  versus  G  ••  o6=l5° 
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FIGURE  6 .31  -  II 

O' jj  r2n/juQn  versus  0  *  o6-30° 
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FIGURE  6.21 -12 

C/y  r2n/juQn  versus  0  :  od  =  45° 
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FIGURE  6. IE  -  13 

C/y  r2n/jjQn  versus  0  :  oL  =  60° 
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FIGURE  6.  EE -I 

m"(o)  versus  n 
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m  {§)  VERSUS  §  :  oi  =  15° 
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FIGURE  6.HT-4 


m  ( (J)  )  VERSUS  (j)  ;  od=45° 
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FIGURE  6.H-5 
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Tr(j)  (o' )  r3n/>u  Qn  versus  n 
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FIGURE  6.  EE-8 

Tr(j)  (o£)  r3n/  /jQn  versus 
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FIGURE  6. HE-  10 
O' ij  r  "n/jjQn  VERSUS  (J)  : 
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FIGURE  6.1H-II 

Oy  r3n/juQn  VERSUS  <j>  J  od  =  30° 
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FIGURE  6. HI-12 

C/y  r3n/;jQn  versus  <j>  ••  o^  =  45° 
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FIGURE  6. HI- 13 

C/jj  r  3n/ju  Qn  versus  §  :  o^=60° 
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FIGURE  6.11-14 


m  (<j))  VERSUS  <j>  :  oC  =  5° 
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FIGURE  6. IE-15 

Gy  r3n/juQn  VERSUS  <j>  :  oC-5° 
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6.2  Error  Analysis  of  the  Numerical  Results 

An  analysis  of  the  errors  involved  in  the  numerical  procedure 
used  is  facilitated  by  the  discussion  presented  in  the  Appendix. 

It  is  shown  in  the  Appendix  that,  if  the  solutions  y^  are  known 
at  x  =  X,  the  accuracy  of  the  step-by-step  integration  procedure  for  n 
first  order  ordinary  differential  equations  is  controlled  by  means  of  a 
test  value  6; 


6 


1  n 

T5  £ 

ID  i  =  l 


(6.1) 


where  y.^  are  the  function  values  at  x  =  X  +  2h,  which  are  obtained 
J  l 

(2) 

by  using  the  increment  Ax  =  h  twice,  and  y^  are  the  function  values 
at  x  =  X  +  2h  as  calculated  in  one  step  with  the  increment  Ax  =  2h.  The 
coefficients  a^  are  the  error  weights,  which  have  to  be  chosen  such  that 


n 

l  a.  =  1  . 

1=1 

The  equations  governing  the  pseudo-plastic  flow  problems  were 
shown  to  consist  of  three  first  order  ordinary  differential  equations,  both 
for  the  plane  and  the  axially  symmetric  flow.  In  order  that  the  local 
truncation  errors  in  the  function  values  y.  are  of  equal  weight,  the 
error  weights  were  chosen  to  be 


'  .  . ■  ■  *  *«  ■  -  -  '  v.  isbi o  .  't 
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od  oJ  froiocb  siaw  aJrigrsw  10*119 


The  maximum  truncation  error,  therefore,  which  could  possibly  occur  in 
the  required  function  at  x  =  X  +  2h,  as  obtained  in  two  steps  with  an 
increment  Ax  =  h  for  each  step,  becomes  from  equations  (A. 25)  and  (6.1) 
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2e  =  6. 


-6 

This  test  value  6,  for  the  control  of  accuracy,  was  chosen  to  be  10 

for  the  execution  of  the  program.  The  maximum  truncation  error  at  each 

-6 

station  in  the  sequence  of  solution  is  therefore  10/2;  consequently, 
the  total  maximum  error  over  the  entire  range,  that  is  from  <J>  =  0  to 
<(>  =  a,  is 


r  .  N  x  10'6 

Total  2 


(6.2) 


where  N  is  the  number  of  stations  needed  to  reach  <J>  =  a.  As  an  example, 
the  solution  to  the  axially  symmetric  converging  flow  problem  with 
a  =  30°,  and  n  =  2/3,  required  128  stations,  hence  substitution  of 


N  =  128 


in  equation  (6.2)  yields  for  the  total  maximum  error  in  the  required 
function,  for  this  particular  case, 

EIotal  ■  6‘4  x  10'5 


no  J6DHU  )  im  Xfifl  f  ir  .ms-',  onq  ertt  nottuo9X9  ariJ  -fo^ 
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over  the  entire  range.  The  errors  involved  in  the  numerical  procedure 
indicate  how  closely  the  exact  solution  is  approximated;  however,  no 
matter  how  good  this  approximation  may  be,  the  errors  do  not  indicate 
the  validity  of  the  solution  in  a  physical  sense.  The  validity  of  a 
solution,  as  applied  to  a  physical  problem,  is  dictated  by  the  mathe¬ 
matical  model  used,  and  the  approximations  made  in  order  to  obtain  that 
solution. 


6.3  Objections  to  the  Power  Law  Model 

The  behavior  of  the  pseudo-plastic  material,  which  is  the 
mathematical  model  considered  in  this  thesis,  is  based  on  the  generali¬ 
zation  of  the  Ostwald's  power  law  model,  and  is  given  by  equation  (1.6), 
which  is 


s 


•  • 
U 


n-1 


(1.6) 


As  mentioned  previously  in  this  thesis,  Reiner  (24)  has  raised  major 
objections  to  the  use  of  this  model.  The  objections,  regarding  the  be¬ 
havior  of  the  apparent  coefficient  of  viscosity,  are 

(i)  for  zero  strain  rates,  the  apparent  coefficient  of  viscosity 
is  infinite,  and 

(ii)  for  infinite  strain  rates,  this  coefficient  becomes  zero. 

This  behavior  of  the  coefficient  of  viscosity,  in  its  limiting  cases, 
can  be  shown  by  consideration  of  FIG.  4.1.  In  this  plot,  the  slope  of 
the  flow  curve  is  a  measure  of  the  apparent  coefficient  of  viscosity,  and 
from  equation  (4.10)  there  obtains 


ohsr  jrt  <’  •  ■•9v  ovr  '  <  :,h  HT  *i  £.-<r  •,->  >n.1  -tsvo 
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n-1 


where,  as  previously  mentioned,  I  and  J  are  the  strain  rate  and  stress 
invariants,  respecti vely .  The  case  n  =  1,  corresponding  to  the  Newtonian 
viscous  liquid,  reduces  to 


dj  _ 
dl  u’ 

where  y  is  the  Newtonian  coefficient  of  viscosity.  However,  for  values 
of  n  in  the  range 


0  <  n  <  1  , 

this  equation  gives  an  infinite  coefficient  of  viscosity  for  1=0,  and 
for  I  =  oo  this  coefficient  becomes  zero.  Therefore,  it  must  be  concluded 
that  for  application  of  the  solutions  to  any  physical  problem,  the  power 
law  model  is  useful  only  for  a  certain  range  of  I;  this  range  must  be 
chosen  such  that  the  limiting  cases,  and  thus  these  objections,  are  ex¬ 
cluded. 

In  the  application  of  the  axially  symmetric  solutions  to  wire 
drawing  processes,  the  objection  of  an  infinite  viscosity  is  removed  by 
the  assumption  that  the  material  does  not  yield  until  some  critical  point 
on  the  flow  curve  is  reached;  this  assumption  has  already  been  discussed 
in  CHAPTER  IV.  Furthermore,  since  the  flow  of  materials  in  experimental 
investigations  is  finite,  the  other  objection  is  removed  as  well. 
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The  two  objections  dealt  with  so  far  are  termed  by  Reiner  (24)  as  the 
"zero"  and  "infinity"  objections,  respectively.  As  previously  mentioned 
in  CHAPTER  I,  an  obvious  defect  in  the  power  law  model  is  that  the  dimensions 
of  y  depend  on  the  value  n.  This  defect,  termed  by  Reiner  as  the  "dimension" 
objection,  has  no  serious  draw-backs  if  experiments,  dealing  with  curve 
fitting,  are  carried  out  for  the  particular  material  under  consideration. 

Although  objections  were  raised  to  the  use  of  the  power  law 
model,  this  model  is  nevertheless  of  importance  in  certain  physical  appli¬ 
cations,  such  as  the  extrusion  of  plastics. 

6.4  Wire  Drawing  and  Hydraulic  Extrusion 

In  CHAPTER  IV  the  velocity  and  stress  fields,  obtained  in 
CHAPTER  III,  for  axially  symmetric  converging  flow  in  a  circular  cone 
were  applied  to  wire  drawing  and  hydraulic  extrusion. 

The  application  of  solutions,  based  on  the  power  law  model,  to 
wire  drawing  seems  to  be  unjustifiable,  since  this  model  describes  a 
material  which  flows  under  any  anisotropic  stress,  no  matter  how  small; 
nevertheless,  it  is  of  academic  interest  and  thus  merits  investigation. 
Furthermore,  since  an  assumption  is  made,  stipulating  a  "yield  limit" 
for  this  pseudo-plastic  material,  thus  rendering  the  so  obtained  model 
as  a  rate  dependent  material  with  a  definite  yield  criterion,  the  results 
of  wire  drawing  become  of  physical  interest.  These  results  become  even 
more  important  if  experimental  investigations,  involving  the  fitting  of 
flow  curves  for  rate  dependent  materials,  are  carried  out.  The  appli¬ 
cation  to  hydraulic  extrusion  does  not  meet  with  this  "yielding"  problem, 
and  should  be  of  major  importance  in  the  extrusion  of  plastics. 
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The  assumption  was  made  that  yielding  does  not  occur  until  the 
stress  invariant  J  reaches  some  value  J  ,  which  is  point  P  on  the  flow 
curve,  see  FIG.  4.2.  The  part  of  the  flow  curve  from  the  origin,  0,  to 
this  critical  point  P,  with  coordinates  IQ  and  J  /y,  could  therefore 
be  considered  as  the  behavior  of  the  material  before  plastic  flow  of 
the  entire  region  ABCD,  see  FIG.  4.3,  first  becomes  possible.  Since 
in  this  thesis  these  preliminary  effects  are  not  taken  into  consideration, 
it  is  immaterial  how  the  flow  curve  behaves  until  that  critical  point  P 
is  reached,  and  hence  the  assumed  flow  curve,  as  shown  in  FIG.  4.2,  should 
be  acceptable. 

Since  the  experimental  investigation,  dealing  with  hydraulic 
extrusions,  was  carried  out  with  a  die  with  a  semi -angle  of  5°,  the 
velocity  profile  and  stress  distributions  for  the  axially  symmetric 
converging  flow  for  that  semi-angle,  with  n  =  0.1,  are  presented  in 
FIGS.  6 . 1 1 1 - 1 4  and  6 . 1 1 1 - 1 5  respectively. 

In  the  application  of  the  pseudo-plastic  flow  fields  to  wire 
drawing,  the  maximum  reductions  are  dictated  by  the  condition  that  yield¬ 
ing  does  not  occur  in  the  drawn  material.  The  results  of  the  relation¬ 
ships  between  the  so  obtained  maximum  reductions  and  semi -angle  a  are 
shown  in  FIG.  6.A-1.  This  investigation  of  the  maximum  reductions  alone 
may  not  be  conclusive  however,  the  fact  that  the  assumption  is  made  that 
the  material  adheres  to  the  walls  of  the  die  may  result  in  negative  normal 
stresses  at  the  die,  if  the  reductions  become  too  large.  The  numerical 
results  show  that  for  any  reduction,  and  z-|  fixed,  the  die  pressure  at  A, 
see  FIG.  6.1,  decreases  as  zQ  increases;  the  die  pressure  remains  positive 
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there  until  =  Z  is  reached,  for  which  the  die  pressure  at  A  is  zero. 

A  further  increase  in  zQ,  and  hence  in  the  reduction,  causes  the  die 
pressure  at  A  to  change  sign,  thus  resulting  in  a  negative  normal  stress 
at  that  point.  This  condition  is  not  physically  tenable,  and  hence 
another  criterion  for  the  maximum  reduction  is  obtained. 


z 


FIGURE  6.1 

Possible  Negative  Normal  Stresses  at  the  Walls 

The  relation  for  the  maximum  reduction,  based  on  the  condition  that 


die 


=  0, 


z=z. 
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is  facilitated  by  equation  (4.15),  which,  using  equations  (4.14)  and 
(4.20)  gives 


2_ 

( 2-3n  1  (nm 1  1  <;i  rub  +  m'rn^rh]-  3n 


as  the  maximum  reduction  which  ensures  that  no  negative  normal  stresses 
occur  at  the  die.  The  results  of  these  maximum  reductions  are  presented 
for  various  semi-angles  in  FIG.  6.A-1.  Clearly,  from  the  numerical  re¬ 
sults  presented  in  FIG.  6.A-1,  the  maximum  reduction  based  on  the  condition 
that  no  negative  normal  stresses  occur  at  the  die  dictates  the  maximum 
reductions  possible.  It  is  interesting  to  note  that  the  special  case 
n  =  0,  considered  by  Shield  (31),  gives  a  maximum  reduction  of  63%  for 
any  semi-angle,  if  the  shear  stress  is  identically  zero  throughout 
the  entire  flow  region;  this  maximum  reduction  is  based  on  the  condition 
that  yielding  does  not  occur  in  the  drawn  material.  It  can  be  shown  that, 
for  this  particular  case,  the  die  pressure  is  positive  along  the  entire 
length  of  the  walls.  However,  Shield  did  not  investigate  the  possibility 
of  negative  normal  stresses  at  the  die,  not  even  for  the  case  where  t 
does  not  vanish  identically  in  the  region  of  flow. 

The  relationship  between  the  ram  pressure  and  semi-angle  a 
for  various  reductions  is  presented  in  FIG.  6.A-3  for  the  hydraulic  ex¬ 
trusion  process.  Similar  results  for  the  wire  drawing  process  have  been 
presented  in  FIG.  6.A-2,  where  in  this  case  the  drawing  stress  is  con¬ 
sidered;  the  maximum  reductions  based  on  the  two  criteria,  discussed  pre¬ 
viously,  are  also  shown  in  this  figure.  It  is  noted  from  FIGS.  6.A-2  and 


. 
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6.A-3  that  both  the  drawing  stress  and  the  ram  pressure  are  monotoni cal ly 
decreasing  functions  of  the  semi-angle  a.  However,  Evans  and  Avitzur  (40) 
found  that  for  a  rigid  perfectly  plastic  solid  there  exists  a  definite 
minimum  for  both  the  drawing  stress  and  the  ram  pressure  as  functions  of 
the  semi -angle  a.  This  minimum  is  denoted  by  the  optimum  semi -die 
angle,  and  is  that  semi -angle  for  which  both  the  drawing  stress  and  the 
ram  pressure  are  a  minimum.  This  distinct  difference  can  probably  be 
attributed  to  the  fact  that  the  material  under  consideration  in  this 
thesis  is  rate  dependent  to  a  significant  extent.  It  would  be  of  interest 
to  investigate  the  relationship  between  the  ram  pressure  and  semi-angle  a 
for  extrusion  of  plastics,  since  these  materials  are  rate  dependent. 


. 

anorJomtf  26  9Tu229‘iq  rnfi'i  9f  J  bni  229*112  pnrwB'fb  art*  rllod  *tol  muml’nrm 

. 

.Insfxs  ins'jr  frprz  6  ol  Jns  >n9q9b  916*1  z\  2  asrtl 

l  .16 - rfTl92  bn6  9'1U229^q  (T16'l  9fil  fl99Wt9<J  q<^  ..oiic'e'l  9f1^  9^6Bf^29Vn  ol 


RESULTS 


OF 

WIRE  DRAWING 
AND 


HYDRAULIC  EXTRUSION 


100 

75 

50 

25 

0 


138 


n  =  o.i 


MAXIMUM  REDUCTION  BASED  ON 
THE  CONDITION  THAT  NO  YIELDING 
OCCURS  IN  THE  DRAWN  MATERIAL 


FIGURE  6  .  A  -  I 

MAXIMUM  PERCENTAGE  REDUCTION  versus  o L 
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6 . 5  Suggestions  for  Future  Work 

Various  results  have  been  presented  for  the  pseudo-plastic 
flow  problems,  considered  in  this  thesis,  however  a  solution  for  the 
Bingham  solid  was  not  obtained.  It  is  of  importance,  therefore,  to 
discuss  a  method  which  could  possibly  yield  an  approximate  solution  to 
the  flow  of  a  Bingham  solid  in  a  circular  cone,  To  this  end  this  chapter 
is  concluded  with  a  brief  discussion  of  this  method. 

It  was  shown  in  CHAPTER  II  that  no  radial,  steady,  quasi-static 
flow  solutions  exist  for  the  plane  flow  of  a  Bingham  solid,  without  body 
forces,  in  a  converging  channel,  or  for  the  axially  symmetric  converging 
flow  in  a  circular  cone,  at  least  if  the  solid  is  assumed  incompressible. 

The  limiting  cases  of  the  Bingham  solid,  however,  the  Newtonian 
viscous  liquid  and  the  perfectly  plastic  von  Mises  solid,  do  indeed  yield 
radial  flow  solutions  to  the  problems  mentioned  above.  Furthermore,  it 
is  clear  that  in  a  long  channel  or  cone  the  velocity  profile  near  the 
entry  section  should  closely  approximate  that  of  the  perfectly  plastic 
von  Mises  solid,  whereas  close  to  the  exit  section,  the  viscous  terms 
dominate  and  hence  the  solution  there  should  be  approximately  that  of  the 
Newtonian  viscous  liquid. 

It  is  the  purpose  of  this  concluding  section  to  outline  briefly 
a  method  which  might  yield  an  approximate  solution  for  the  flow  of  a 
Bingham  solid  in  a  circular  cone  of  small  semi-angle  and  of  finite  length. 
This  method  was  successfully  applied  by  Sutterby  (37)  to  the  axially 
symmetric  converging,  steady,  incompressible  Newtonian  viscous  flow  in  a 
circular  cone  of  finite  length,  and  of  small  semi-angle,  with  specified 
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conditions  at  the  inlet  section. 

Let  the  flow  of  a  Bingham  solid  in  a  circular  cone,  of  length 
L,  be  described  by  spherical  polar  coordinates,  (r,<j>#6);  the  z-axis  is 
taken  along  the  axis  of  symmetry,  see  FIG.  6.2.  Further  necessary 
notation  is  also  shown  in  this  figure. 


FIGURE  6.2 

The  Cone  of  Finite  Length 


Since  no  radial  flow  exist  for  the  Bingham  solid,  the  flow  must 
consist  of  a  tangential  component  as  well.  The  flow  field,  therefore, 
consists  of  a  radial  velocity  vr  =  v  ( r ,4>) ,  a  tangential  velocity 

v,  =  v.(r,<|>),  and  vn  =  0.  From  equations  (2.1),  there  obtains  for  the 

(p  (p  u 
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components  of  the  strain  rate  tensor 


3v. 


j  _  _ r_ 

r  ~  3r  » 


d  + 

<j>  r  r  34>  * 


d„  = 


v,cot<f> 


0  r 


+  :* 


(6.3) 


T  3v  *j  3v , 

d  =  1 _ l  +  1  _i 

r4>  2r  3<t>  2  Sr  ' 


The  equation  of  continuity  becomes 


13  /  y  2  \  1  3  /  *j.\ 

?3?(r  vr  +  Hri*  9*  (Vin*) 


=  0 


(6.4) 


The  strain  rate  invariant, 


I 


/u 


mn 


(6.5) 


is  obtained  from  the  strain  rate  relations  (6.3).  Similarly,  the  stress 
deviator  components  are  found  by  substitution  of  equations  (6,3)  and  (6,5) 
in  the  constitutive  equations  (1,4).  These  expressions  are  very  complex, 
and  will  not  be  reproduced  here.  The  equations  of  motion,  assuming  steady 
flow,  and  without  body  forces,  are 


Sn  9V 
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9v ,  9v,  9s,  9s  , 

iE.+  D(rv  _±+v  _i)=— i.+  r— H 
3<j)  pv  vr  3r  v<{.  34>  J  3<J>  r  3r 


+  3s 


r<j) 


(6.6) 


Substitution  of  the  stress  deviators  in  equations  (6.6),  results  in  a 
set  of  equations,  which  along  with  equation  (6.4)  constitute  three  equations 

for  the  three  unknowns  p,  v  .  and  v,  . 

r  <J> 

If  the  circular  cone  under  consideration  is  of  small  semi-angle 
a,  then  it  may  be  assumed  that  the  streamlines  are  only  slightly  curved. 
Consequently,  v^  is  small  compared  to  v^,  and  furthermore,  sin<j)  and  cos<t> 
may  be  replaced  by  4>  and  1,  respectively.  With  these  approximations,  the 
governing  equations  are  simplified  to  some  extent. 

The  resulting  equations  of  motion  are  of  the  following  form 


fjr  +  pF  +  kG  +  pH  =  0  , 


||  +  pF  +  kG  +  yH  =  0  , 


and  the  equation  of  continuity  becomes 


9?”  +  2vr  +  a/  +  /  =  0  ’ 


(6.7) 


(6.8) 


where  the  functions  F,  G,  H,  7,  G,  and  FT  are  functions  of  r  and  <t>,  and 
of  the  velocity  components  and  their  derivatives  with  respect  to  r  and  <j> 
up  to  and  including  their  second  order.  The  volume  flow  condition,  for 
small  semi -angles,  becomes 
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[CL 

vr  (j)d(j)„ 

< 

0 


(6.9) 


The  problem  reduces  to  the  solution  of  equations  (6.7)  and 
(6.8)  subject  to  the  volume  flow  condition  (6.9),  and  the  following 
boundary  conditions: 
from  symmetry, 


9cj)  vcj)  "  °»  at  cj)  -  0, 

and  since  there  is  no  flow  across  the  walls  of  the  die, 

v^  =  0  at  <J)  =  a, 

and  finally  the  assumed  no-slip  boundary  condition 

-  0,  at  <f>  =  a. 

In  this  method  the  velocity  profile  at  the  entry  section,  that 
is  r  =  R  ,  is  assumed  to  be  that  of  the  Bingham  solid  in  a  circular  tube, 
and  the  pressure  is  assumed  to  be  constant  there.  Consequently,  the  in¬ 
let  conditions  are 

at  r  =  Rq:  v^  =  0,  p  =  pQ  ,  and  vr  is  as  given  by  the 


Bingham  solution  in  a  circular  tube.  A  finite  difference  analysis  can 


, 
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now  be  applied  to  the  region  of  deformation. 

It  is  noted  that  the  complexity  of  the  system  is  more  pronounced 
than  that  for  the  Newtonian  viscous  flow.  This  may  be  a  defect  in  this 
method,  as  applied  to  flow  of  materials  with  a  less  simple  constitutive 
equation.  A  more  thorough  theoretical  investigation  is  therefore  re¬ 
quired;  if  a  solution  is  obtained,  these  results  should  be  of  major 
importance  for  hydraulic  extrusion  and  wire  drawing  processes. 
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APPENDIX 


NUMERICAL  INTEGRATION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 
A. 1  Introduction 

The  governing  equations  of  CHAPTERS  III  and  IV  were  integrated 
by  Gill's  variation  of  the  Runge-Kutta  fourth  order  method  (41).  This 
is  a  self  starting  integration  procedure  which  is  applicable  to  a  system 
of  n  first  order  ordinary  differential  equations.  The  following  expla¬ 
nation  first  considers  a  single  differential  equation,  this  is  then 
generalized  to  a  system  of  ordinary  differential  equations,  and  finally 
the  accuracy  of  the  method  is  indicated. 

A. 2  The  First  Order  Ordinary  Differential  Equation 

Consider  the  single  first  order  ordinary  differential  equation 


The  initial  condition  is  y  =  Y  at  x  =  X  and  suppose  that  y  is  required 
at  x  =  X  +  h,  where  h  is  small.  One  possible  procedure  is  to  consider 
the  Taylor  series  expansion  about  the  point  P,  with  coordinates  (X,Y), 


y  -  Y  =  6y  -  h($4  +  jj-  (^4)  +  ft,-  (^-4) 


dx'p  21  Vdx2,p  31  'dx3  P 


(A.  1 ) 


where  the  coefficients  are  calculated  from 
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&  = 
dx 


f(x,y) 


9 


d2.y  -  3f  I  f  3f 

dx2  3x  3* 


etc. 


The  complexity  of  the  coefficients  renders  this  process  impracticable. 

An  alternative  procedure,  introduced  by  Runge  and  modified  by 
Kutta,  proceeds  as  follows: 

Define,  following  Gill's  notation, 


k0  =  hf (X,Y)  , 
k-j  =  hf ( X  +  mh,  Y  +  mkQ)  , 
k^  =  hf ( X  +  nh,  Y  +  (n-r)kQ  +  rk-j)  , 
k^  =  hf ( X  +  ph,  Y  +  (p-s-t)kQ  +  s k-j  +  tk^)  , 

and  construct 


6y  =  akQ  +  bk-j  +  ck^  +  dk^  .  (A. 2) 

FIGURE  A.l  indicates  the  geometrical  nature  of  this  procedure.  The  six 

constants  m,  n,  p,  r,  s,  t,  and  the  four  weighting  factors  a,  b,  c,  and  d 

/  4\ 

are  chosen  such  that  the  terms  up  to  and  including  0(h  )  in  the  Taylor 
series  expansion  of  equation  (A.l)  agree  with  the  corresponding  terms  in 
the  expansion  of  equation  (A. 2). 


,9fd6:  rioi  Kiur  aasooiq  ai.riJ  a<r‘  i^n  aim  bi  l  ts.-  sr  \  M  ^  oo  eri 

- 

)efttboiT'  bne  sorwH  yd  uSDUbcn:  i  ,  9iub9jO',q  ov  *MSr  t-  <•<; 
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FIGURE  A. 1 

Geometrical  Construction  of  6y 


A. 3  Equations  Governing  the  Constants  and  Weighting  Factors 

The  Taylor  series  expansion  of  equation  (A. 2)  requires  the  ex¬ 
pansions  of  the  weighting  functions  kQ,  k-j ,  k^,  and  k^.  These  expansions 
are  obtained  by  making  use  of  the  Taylor  series  expansion  in  two  variables 
(42) 


—2  —3 

f(x+p,y+q)  =  f(x,y)  +  Df(x,y)  +  jy  f(x»y)  +  jv  f(x»y)  + 
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where  the  operator  D  is  given  by 

Br-<P!x  +  0fy>r- 

With  the  introduction  of  the  notation 

f0  =  f(X,Y), 

the  expression  for  k0  becomes 


k  =  hf  . 
o  o 


Letting 


D  =  —  +  f  — 

3x  o  3y  * 


then 


mtl  h  +  mkO  ly  =  mh  <lx  +  fo  fy* 


consequently 


k,  = 


2.  2n2,  3u3n3, 

,  rf  .  ,  ,  m  h  D  f  ,  m  h  D  f 

h  [f  +  mhDf  +  — 2l - +  — 31 - H 


To  evaluate  k^,  then 


=  mh  D, 


...] 

p 


nh  lx  +  Kn-r>k0  +  rkl]  ly  =  nhD  +  r(krkO)  ly  = 


nor"  'ton  s.tt  to  nor  Jo  jo  1  r  /rV; 
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nhD  +  mrh2  [Df  +  +  m-Vp3f  +  . . .  ]  |_  , 

c  .  o .  p  9y 


consequently 


2  2  2  3  3  3 

n  n  D^f  ,  nJhVf 


k2  =  h[f  +  nhDf  +  11  1  + 


3! 


+  ... 


+  mrh2  {f  Df  +  5y  f  D2f  +  nhDf  •  Df  +...}]  , 

y  y  y  p 


where 


f  -Si 

y  ay 


Similarly  to  evaluate 


ph  37+  [(p-s-t)kQ+sk1+tk2]  =  phD  +  [s(k1  -kQ)  +  t(k2-kQ)]  |p-  = 


phD  +  [h  ms  {Df  +  — ^T”  + 


2  2  2  3 

mhD  t  ,  mTDJf 


3! 


+  . . . } 


+  h2nt  {Df  +  nhD^f  +  nWf  +  _} 


2! 


3! 


+  h3tmr  (f  Df  +  S-  f  D2f  +  nhDf  •  Df  +...}]  |-  , 


2!  y 


3y 


consequently 


2h2n2f  3,3  « 

k,  =  h[f  +  phDf  +  +  Syf-  DJf  +  ... 


3! 


+  h2(msDf+tnDf)f  +  jy  (sm2D2f  +  2tmrf  Df  +  tn2D2f)  f 


+  to  •  tGrln  +  tS0  1  +  to  t}  SfHm  + 


...  +  V'd  +■  +  tGriq  +  t]ri  «  i 
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+  n  (psmDf  +  ptnDf)  Df  +  ...]  . 


The  expansion  of  equation  (A.l)  becomes 


2  3 

6y  =  [hf  +  {q-  Df  +  yy  (D2f  +  f  Df)  + 


xr  (D3f  +  f  D2f  +  f  2Df  +  3DfDf  )  +  . 
4!  y  y  y 


x=X 


(A. 3) 


The  equations  connecting  the  constants  and  the  weighting 
factors  are  obtained  by  equating  the  eight  terms  in  the  expansion  of  6y, 
equation  (A. 3),  to  the  corresponding  terms  in  the  expansion  of  6y  =  akQ  + 
bk i  +  ck^  +  dk^  .  There  results  the  equations 


a  +  b  +  c  +  d=  l  , 


bm  +  cn  +  dp  =  1/2  , 
bn/  +  cn^  +  dp^  =  1/3  , 

(A. 4) 

bn/  +  cn^  +  dp^  =  1/4  , 
crm  +  dsm  +  dtn  =  1/6  , 
crn/  +  dsn/  +  dtn^  =  1/12  , 


drmn  +  dsmp  +  dtnp  =  1/8  , 


dtrm  =  1/24  . 


The  set  of  equations  (A. 4)  represent  eight  constraints  on  the 
ten  unknowns  a,  b,  . ..,  s,t,  consequently  two  further  consistent  relations 
may  be  imposed;  this  may  be  done  in  various  ways,  but  first  it  is  shown 
that  p  =  1 . 

A. 4  Determination  of  the  Constants  and  Weighting  Factors 

Adding  the  second  of  the  set  of  equations  (A. 4)  multiplied  by 
mp  and  the  third  multiplied  by  -  (m  +  p)  to  the  fourth  equation,  gives 


cn(m-n)(p-n)  =  f..  +  l  . 


(A. 5) 


Multiplying  the  fifth  equation  by  p  and  subtracting  from  this  the  seventh 
equati on , 


cmr(p-n)  =  5  -  ^  • 


(A. 6) 


Similarly  from  the  fifth  and  sixth, 


dnt(n-m)  =  jg  "  ^  • 


(A. 7) 


Eliminating  d  between  equation  (A. 7)  and  the  eighth  of  the  set  of  equations 


. 
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Substituting  equation  (A. 8)  in  equation  (A. 6), 


cn(p-n) (m-n)  =  (2m-l)(|  -  |) 


(A. 9) 


Comparing  equation  (A. 9)  with  equation  (A. 5)  gives 


mp  =  m  . 


(A. 10) 


From  the  eighth  of  the  set  of  equations  (A. 4) 


m  f  0  , 


hence  the  solution  to  equation  (A. 10)  reduces  to 


P  =  1 


The  values  of  a,  b,  c,  and  d  are  determined  from  the  first 
four  of  the  set  of  equations  (A. 4), 

a  +  b  +  c  +  d=  l  , 


bm  +  cn  +  d  =  j  , 


2  2  1 

bm  +  cn  +  d  =  y  , 


(A. 11) 


bm^  +  cn^  +  d  =  , 


r 

.  -  q 


since  p  =  1.  This  set  of  equations  yields  a  unique  solution  for  a,  b, 
c,  and  d,  if  the  coefficient  determinant  is  non-zero  that  is  if 
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mn(m-n)(n-l)(l-m)  1  0  . 


(A. 12) 


Solving  the  set  of  equations  (A. 11); 


a 


1  j.  1  -  2(m+n) 

2  12mn 


k  -  2n-l 

T2m(n-m)(l-m)  * 

(A. 13) 

.  _  1 -2m 

c  12n(n-m) ( 1-n)  5 


1  x  2(m+n)  -  3 

2  12(T-m) (1-n) 


Finally,  the  constants  r,  s,  and  t  are  determined  from  the 
fifth,  sixth  and  seventh  of  the  set  of  equations  (A. 4), 

crm  +  dsm  +  dtn  =  , 

crn/  +  dsm^  +  dtn^  =  »  (A.  14) 

crmn  +  dsm  +  dtn  =  . 


Once  again  a  unique  solution  is  obtained  if  the  coefficient  determinant 
of  the  set  of  equations  (A. 14) 


* 


■■ 


J:\ib  +  rite  ♦  no 
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cd2m2n(n-m)(n-l)  f  0  .  (A. 15) 

Consequently,  the  unique  solution  to  r,  s,  and  t  is 

r  -  n(n-m) 
r  '  2m(l-2m)  » 

=  (l-m){m+n-l-(2n-l )2} 

2m(n-m) {6mn-4(m+n)+3}  * 

t  =  (l-2m)(1-m)(l-n) 

n(n-m) {6mn-4(m+n)+3}  * 

The  parameters  m  and  n  may  be  assigned  in  any  manner  consistent 
with  the  previous  equations.  Let  the  interval  h  be  divided  into  three 
equal  parts,  then 


(A. 16) 


Substitution  of  equations  (A. 16)  in  the  set  of  equations  (A. 13)  yields 


Similarly  the  constants  r,  s,  and  t  are 


r  =  1  , 

s  =  -  1  , 

t  =  1  . 


The  weighting  functions  obtained  are 


kQ  =  hf (X,Y) , 

h  kn 

k1  =  hf (X  +  f  •  V  +  3^)  , 

k2  =  hf(X  +  fl.Y  -  ki)  , 

k,  =  hf (X  +  h,  Y  +  k  -  k,  +  k,)  , 
d  o  1  2 


and  the  increment  in  y  becomes 


6y  =  l  (k  +  3k] 


+  3)<2  +  k^)  . 


The  formulae  for  the  k's  and  the  increment  in  y  thus  obtained  are  those 


.(Y.XVM  *  .  * 


(5  +  Y  ,  +  X)Yrt  »  * 
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due  to  Kutta  (42) . 

Another  useful  choice  for  the  parameters  m  and  n  is  facilitated 
by  setting  the  determinants  (A. 12)  and  (A. 15)  equal  to  zero.  The  case 
of  interest 


m  =  n  (A. 17) 

is  considered.  Substitution  of  equation  (A. 17)  in  equations  (A. 13) 
further  imposes  the  conditions 


in  order  that  b  and  c  remain  finite.  Substituting  these  values  for  m 
and  n  into  the  first  and  fourth  of  the  set  of  equations  (A.ll)  yields 


Either  b  or  c  is  now  arbitrary  but  not  both,  since  the  condition 


a  +  b  +  c  +  d=  l 


r 34  .  •  \i  (V  •)  n  i  tc  n:  tul  ;;  u i  ..  e^-^  'rzr.o:  ?\ 


•'  M  Ogin  'v*rlj  , 


nr>  nfBW  o  bn  a  d  dfirid  tabio  n'r 


is  to  be  satisfied.  Assuming  that 


equations  (A. 14)  give 


s  =  0  , 

t  =  1  . 

The  constants  and  weighting  factors  thus  obtained  give  rise  to  a  set  of 
formulae,  originally  due  to  Runge  (42),  for  the  k's  and  an  increment  in 

kQ  =  hf(X.Y)  , 

h  kn 

k,  =  hf(X  +  \  ,  V  +  f-  )  , 

h  kl 

k,  =  hf (X  +  Y  +  j-)  , 


k,  =  hf (X  +  h,  Y  +  k2)  , 


,  (  ji  +  V  ,rf  +  X)t.1  =  i 
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and  the  increment  in  y  becomes 

6y  =  g-  (k0  +  2k1  +  2k2  +  k3)  . 

Two  methods  of  evaluation  of  the  constants  have  been  considered, 
namely  those  due  to  Kutta  and  Runge.  Another  approach,  due  to  Gill  (41), 
is  based  on  the  number  of  storage  registers  needed  in  the  calculation  of 
each  step.  Gill's  variation  will  now  be  considered  by  application  to  a 
system  of  ordinary  differential  equations. 

A. 5  The  System  of  First  Order  Ordinary  Differential  Equations 

The  consideration  of  the  system  of  equations  is  facilitated  by 
the  preceeding  discussion  of  a  single  first  order  ordinary  differential 
equation.  The  system  of  equations  is 

dyi 

d3T  =  VX,V  ;  1  ,J  =  1  j2,  •" ,n  * 

with  the  initial  values  of  Y,,  Y2,...,Y  at  x  =  x-  The  weighting  functions 
for  the  system  are  expressed  as 

kio  =  hfi(X,Yl,Y2’  ’ 

k^i  =  h  f  ^  ( X  ^  mh ,  V-|  mk^Qj  Y2  mk^Q*  •••)  * 
k-2  =  hf. (X  +  nh,  Y,  +  (n-r)k1Q  +  rk^,  Y£  +  (n-r)k2Q  +  rk21 ,  ...)  , 


.be.abtanoo  nead  avert  adnedanoo  artd  *o  nor*eufeva  *bortJan,  ow  'J| 
\0  norJefuofeD  arid  nt  babaan  aiadar^  sBeioda  *>  i*mm  *11  no  baaed  a' 
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ki3  =  h1VX  +  ph>  Y1  +  (P-s-t)k-jo  +  skn  +  tk22>  ...)  , 

and  the  fourth  order  approximation  is 

6yi  =  akio  +  bki 1  +  cki2  +  dki3  * 

Gill  has  introduced  a  notation  which  simplifies  the  theory  and  is  also 
better  adapted  to  computer  programming  techniques.  Let 

x  =  yQ  ,  (A. 18) 

then  the  system  of  equations  becomes 

dy  • 

-  f  .j  ( y  j )  »  i  ~  1 » 2 »  . .  • ,  n , 

J  —  0>1»  •  •  <  »n  j 

where  f  =1.  Under  the  transformation  (A. 18)  the  weighting  functions 
o 


k 


i2 


kil  =  hfi(Yj  +  mkjo}  > 

=  hfi(Yj.  +  (n-r)kj0  +  rk.,) 


9 


are  obtained  as 


' 


.  (tnt«  *  o : 


l  vV.M  =  i 
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ki 3  =  hfi(Yj  +  (p-s-t)kj0  +  s k j -j  +  tkj2)  , 

where  i  =  1,2,  ...,  n  ,  and  j  =  0,1,  ...»  n. 

The  Taylor  series  expansion  for 

<5yn-  =  y.j(x  +  h)  -  y i ( x) 

about  the  initial  condition  is  now  obtained.  The  manipulations  are  con¬ 
siderably  simplified  by  the  introduction  of  the  notation 


fi  ■  VV  * 


and 


f.o  =  (!!i) 

1  ^j(X.Yj)  ’ 


(A. 19) 


where 


i  =  1 ,2,  . . . ,  n  , 


j  0,1,  . . . ,  n  . 


The  Taylor  series  expansion,  as  given  by 


<5yi  =  y^X+h) 


dyi 

-yi(x)  =  h(^) 


,2  d2y. 

+  — TT* 

(X,Y.)  ’  dx  (X,Y. ) 


+  n!  ft  . 

31  j,3 

dx  (X,Y.) 


, 


bn6 
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becomes  with  the  notation  of  (A. 19) 


«yi  =  y^x+h)  -  y.(x)  =  hf .  +  y. 


IT  (fjVijk  +  Vi  V> + 


(A. 20) 


L/f  f  f  f  j‘kl  +  3f  f  f  ^  +  f  f  f  ^  +  f  f  Jf  kf  1\  + 

4!  ujTkTlTi  J  jVl  Ti  TjTkTl  Ti  TjTk  T1  Ti  •  * 


Taylor  series  expansions  are  also  required  for  the  weighting 
functions  k^,  k^,  k.g»  and  k^>  the  expression  for  k^Q  simply  reduces 
to 


kio  =  hfi  • 


(A. 21) 


The  expansion  of  k^  is  obtained  from 


df. 


2  2 
rn  k .  k  d  t  . 


fi(Yj  +  mkJo> 


W  +  mkjo(dy|)Y  +  21°  P°  (dyjdyp)y  +  •"  ' 

J  <1 

(A. 22) 


Substituting  equation  (A. 21)  in  the  expansion  (A. 22)  gives 


ki  1  =  hqtYj  +  mkj0)  ■  hfi  +  mh2fjfiJ  +  V"  fjVl 


jk 


+  infill  f  f  f  f  Jkl  + 

6  j  k  1  i 
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The  substitutions  and  manipulations  required  to  find  the  expansions  for 
k.j 2  and  k. ^  are  complicated,  and  it  serves  no  purpose  to  reproduce  them 
here.  However,  it  can  be  shown  that  the  equality  of  the  eight  terms  in 
equation  (A. 20)  and  the  corresponding  ones  in  the  expansion  of  = 
ak.jo  +  bk^.-j  +  ck^  +  dk^  leads  to  the  set  of  equations  (A. 4)  previously 
given  for  a  single  equation.  Table  A.l,  extracted  from  Gill's  paper  (41), 
indicates  the  typical  terms  as  well  as  their  coefficients. 


TABLE  A.l 

Typical  Terms  and  their  Coefficients 


Term 

Coefficient  in 

<$yi=yi(x+h)-yi(x) 

6y  j  ~  ak^.  ^+bk^  -j+ck^  2"^"^^i  3 

hfi 

h2fjV  . 

h3f .fkfijk 

1 

a  +  b  +  c  +  d 

1/2 

bm  +  cn  +  dp 

1/6 

l/2(bm2  +  cn2  +  dp2) 

hVkjfik 

1/6 

crm  +  d(sm  +  tn) 

h4f 

1/24 

l/6(bm3  +  cn3  +  dp3) 

h^f  f  f  ^f  ^ 
n  j  kTl  Ti 

1/8 

crm  +  d(sm  +  tn)  p 

i4.p  .  f  jko  1 

h  fjfkfi  fi 

1/24 

l/2{crm2  +  d(sm2  +  tn2)} 

4.  ,  f  1 

h  fjfkf1  fi 

1/24 

dtrm 

A. 6  Gill's  Variation  of  the  Runge-Kutta  Fourth  Order  Method 

The  procedure  Gill  used  in  removing  the  indeterminancy  of  the 
constants  is  now  discussed;  the  purpose  of  this  procedure,  as  previously 
noted,  is  to  minimize  the  number  of  storage  registers  required  for  auto- 
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matic  machine  computations.  Further,  Gill  chose  the  parameters  such 
that  the  accuracy  was  the  highest  attainable. 

In  starting  out  the  integration  procedure  n  +  1  registers  are 
required  to  store  the  quantities  Y. ,  j  =  0 ,1 ,  . . . ,  n.  At  the  end  of  the 

J 

first  stage  an  extra  register  is  required  to  store  the  quantity  k^Q. 

During  the  second  stage,  the  quantities  Y.  +  mk .  and  k.-,  are  to  be 

0  jo  J 

stored.  Before  proceeding,  however,  the  quantities  Y.  +  (n-r)k.  , 

J  J  ® 

Y.  +  (p-s-t)k.  and  Y.  +  ak.  should  be  stored  as  well  since  they  are 

J  J  ^  J  J 

needed  in  the  following  stages.  These  five  quantities  are  linearly  de¬ 
pendent  and  can  be  represented  by  three.  For  the  third  stage  the  quantities 

Yj  +  (n-r)kjo  +  rkjl’  Yj  +  (P-$-t)kjo  +  skjl  >  Yj  +  akjo  +  bkjl>  and  kj2 
need  to  be  stored,  therefore,  four  registers  are  required.  At  the  final 

stage  the  quanties  to  be  stored  are  three  in  number,  namely 


Y.  +  (p-s-t)k.  +  sk.n  +  tk.0,  Y.  +  ak.  +  bk.,  +  ck.9,  and  k.Q  . 

J  jo  jl  j2  j  jo  jl  j^  jd 

So  in  general  three  storage  registers  are  enough  except  at  the  third 
stage,  where  four  are  needed.  This  draw-back  can  be  overcome  by  making 
the  first  three  linearly  dependent;  this  condition  implies  that 


1  n-r  r 

1  p-s-t  s 

1  a  b 


=  0  , 


which  reduces  to 


■ 


.  ,  +  S  3 


* 
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Solving  the  set  of  equations  (A. 24),  gives 

a  =  1/6  , 

b  -  jo +4i • 

c  ■  |  [i  ±4  3  , 

d  =  1/6  , 

and 

r  =  1  + /i , 
s  =  +  /  j  , 
t  =  1  +  /l  , 

with 

m  =  n  =  1/2  , 

p  =  1  . 

c 

Gill  has  shown  that  in  making  the  error  in  the  term  of  0(h  )  in  the  Taylor 
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[p  -  ( s+t ) ] ( b-r)  +  (n-r)(s-b)  +  a(r-s)  =  0.  (A. 23) 

The  system  of  equations  to  be  solved,  therefore, reduces  to  the  set  of 
equations  (A. 4)  and  in  addition  equation  (A. 23). 

It  has  been  shown  previously  that  p  =  1;  then,  if  m  and  n  are 
constrained  to  be  equal  to  1/2,  a  set  of  equations  is  obtained,  which  for 
completeness  is  shown  below 


a  +  b  +  c  +  d=  l  , 

b  +  c  +  2d  =  1  , 

b  +  c  +  4d  =  4/3  , 

b  +  c  +  8d  =  2  , 

cr  +  ds  +  dt  =  1/3  , 

cr  +  ds  +  dt  =  1/3  , 

cr  +  2ds  +  2dt  =  1/2  , 

dtr  =  1/12  , 

[1  -  ( s+t )] (b-r)  +  (■£--  r)(s-b)  +  a(r-s)  =  0  . 


(A. 24) 


• 

«  d\ f  :  : 


,  X  \  T  *■  z 
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series  expansion  the  least,  and  thus  obtaining  maximum  accuracy,  the 
upper  sign  of  the  square  root  should  be  used.  The  formulae,  due  to  Gill, 
then  become 


kio  h1Vyoo,y10’y20*  *•*)  ’ 
kil  =  ^i  ^y01  ,yl  1  ,y21  *  *  ’  *  ^  * 
k i 2  =  hf.j (^02  ,yl 2 ,y22  *  *  *  *  ^  * 

ki3  =  h1Vy03,y13,y23’  * 


where 


y .  =  Y  •  , 

<yio  i  * 


yil  Yi  +  2  k  ’ 


yi2  =  Yi  +  ?+  ]kio  +  [1  '  ]  kil  > 


y -j 3  =  Yi  +  y  \  ]kil  +  H  +  y  2  ^  ki2  ’ 


and  the  accepted  value  of  y..  at  the  end  of  the  step  is 


yi4  =  Yi  +  F  kio  +  T  [1  '  /  ?  ]kil  +  ^[1  +  /  2  ]ki2  +  6  ki3  ’ 


*  5  +  Sr51,  1  ^  +  r]l  +  U*  1  V  '  ]  I  +  of*  I  +  f  =  M* 
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where 


i  =  0,1,  . . . ,  n  . 

A  subroutine  based  on  Gill's  variation  of  the  Runge-Kutta  fourth  order 
method  (43)  has  been  prepared,  and  is  shown  in  the  body  of  this  Appendix. 

Up  to  now  no  mention  has  been  made  about  the  errors  encountered 
in  using  this  process;  this  is  done  in  the  next  section. 

A. 7  Brief  Discussion  on  the  Accuracy  of  the  Method 

In  order  to  discuss  the  results  qualitatively,  it  is  necessary 
that  some  estimate  is  made  of  the  errors  due  to  rounding-off  and  truncation. 

c 

The  latter  are  caused  by  the  fact  that  terms  of  0(h  )  and  higher  are  ne¬ 
glected  in  the  Taylor  series  expansion  about  the  initial  values. 

Gill  has  discussed  the  rounding-off  errors  of  such  quantities 
as  /  and  ^  at  length.  However,  since  the  program  is  written  in  double 
precision,  that  is  the  calculations  are  performed  carrying  16  significant 
figures,  it  is  reasonable  to  assume  that  the  errors  introduced  due  to 
rounding-off  are  negligible. 

Truncation  errors  or  estimates  thereof  are  not  accounted  for 
in  the  program.  In  order  to  overcome  this  disadvantage,  control  of  accuracy 
is  accomplished  by  adjusting  the  increment  in  x;  a  comparison  is  then  made 
between  the  function  values  obtained  at  x  =  X  +  2h  by  firstly  using  the 
increment  Ax  =  h  in  two  stages  and  secondly  by  using  double  the  increment. 

Ax  =  2h,  in  a  single  stage. 

Returning  now  to  the  single  first  order  equation 


-9n  9-ffi  -i9riprri  bne  (2ri)0  zm$3  dbrtdt  toeT'  srtt  yd  b92UB3  9^6  9riT 

.riigne:  3b  g  bne  26 

.9fdrgrrg9a  9^6  t^o-gnrbnucn 
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&L  = 
dx 


f(x,y) 


Suppose  that  the  initial  condition  is  y  =  Y  at  x  =  X,  and  that  the 
function  value  is  required  at  x  =  X  +  2h,  where  h  is  small.  Let  the 
error,  due  to  truncation,  of  the  function  value  at  x  =  X  +  h  be  denoted 
by  e.  Since  terms  0(h5)  are  neglected  in  the  Taylor  series  expansion, 


e  = 


ch 


5 


approximately,  where  c  is  some  constant.  The  error  in  the  function  value 
at  x  =  X  +  2h,  as  calculated  in  two  stages,  is  then 


Secondly,  the  increment  Ax  =  2h  is  used  and  consequently  the  error  in  the 
function  value  at  x  =  X  +  2h,  as  obtained  in  a  single  stage,  is  given  by 

s'  =  32ch5  , 


consequently 


2e  =  2ch5  =  (e*  -  2e)  .  (A. 25) 


From  the  preceeding  discussion,  this  may  be  written  as 


26  (133Jf-W  9'i  6m  ailIJ  .Kilt  -  >'■  !  3  !  '  f  S'  '  nK>“’ 


2e  ■  ts  [y(2)  -  y(1)  ]  , 
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where  y^  (X  +2h)  is  the  function  value  obtained  in  two  stages  with  the 
increment  Ax  =  h,  and  y^  (X  +  2h)  is  the  function  value  calculated  in 
one  stage  with  the  increment  Ax  =  2h.  This  discussion  of  the  truncation 
errors  may  readily  be  extended  to  a  system  of  ordinary  differential 
equations. 

In  the  subroutine  a  test  value  6  is  generated  for  the  control 
of  accuracy  and  is  defined  as 


6 


1 


n 


T5  I 
IJ>  i  =  l 


5 


where  y.^  and  y.^  have  the  same  meaning  as  before,  and  the  coefficients 
a.  are  the  error-weights.  The  test  value  6  is  an  approximate  measure  of 
the  local  truncation  error  at  the  point  x  =  X  +  2h,  and  may  be  specified 
arbitrarily  at  the  start  of  the  program.  The  test  value  6  is  denoted  by 
DELT  in  the  subroutine.  For  completeness  the  subroutine  based  on  Gill's 
variation  of  the  Runge-Kutta  fourth  order  method  is  now  presented. 


' 

fsttnona^fb  ^i6nhno  meJa^a  6  of  b9bri9fx9  9d  rben  an<ms 
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A. 8  Subroutine  DRKGS 


SUBROUTINE  DRKGS 

PURPOSE 

TO  SOLVE  A  SYSTEM  OF  FIRST  ORDER  ORDINARY  DIFFERENTIAL 

EQUATIONS  WITH  GIVEN  INITIAL  VALUES. 

USAGE 

CALL  DRKGS ( PRMT, Y ,DERY,NDIM,IHLF,FCT,OUTP,AUX) 

PARAMETERS  FCT  AND  OUTP  REQUIRE  AN  EXTERNAL  STATEMENT. 

DESCRIPTION  OF  PARAMETERS 

PRMT  -  DOUBLE  PRECISION  INPUT  AND  OUTPUT  VECTOR  WITH 
DIMENSION  GREATER  THAN  OR  EQUAL  TO  5, WHICH 
SPECIFIES  YHE  PARAMETERS  OF  THE  INTERVAL  AND  OF 
ACCURACY  AND  WHICH  SERVES  FOR  COMMUNICATION  BETWEEN 
OUTPUT  SUBROUTINE  (FURNISHED  BY  THE  USER)  AND 
SUBROUTINE  DRKGS.  EXCEPT  PRMT(5)  THE  COMPONENTS 
ARE  NOT  DESTROYED  BY  SUBROUTINE  DRKGS  AND  THEY  ARE 

PRMT ( 1 )-  LOWER  BOUND  OF  THE  INTERVAL  (INPUT), 

P RMT ( 2 ) -  UPPER  BOUND  OF  THE  INTERVAL  (INPUT), 

PRMT (3) -  INITIAL  INCREMENT  OF  THE  INDEPENDENT  VARIABLE 
(INPUT), 

PRMT (4)-  UPPER  ERROR  BOUND  (INPUT).  IF  ABSOLUTE  ERROR  IS 
GREATER  THAN  PRMT(4)>  INCREMENT  GETS  HALVED. 

IF  INCREMENT  IS  LESS  THAN  PRMT(3)  AND  ABSOLUTE 
ERROR  LESS  THAN  PRMT(4)/50,  INCREMENT  GETS  DOUBLED. 
THE  USER  MAY  CHANGE  PRMT(4),BY  MEANS  OF  HIS  OUTPUT 
SUBROUTINE. 

PRMT (5)-  NO  INPUT  PARAMETER.  SUBROUTINE  DRKGS  INITIALIZES 
PRMT ( 5 ) =0 .  IF  THE  USER  WANTS  TO  TERMINATE 
SUBROUTINE  DRKGS  AT  ANY  OUTPUT  POINT,  HE  HAS  TO 
CHANGE  PRMT (5)  TO  NON-ZERO  BY  MEANS  OF  SUBROUTINE 
OUTP.  FURTHER  COMPONENTS  OF  VECTOR  PRMT  ARE 
FEASIBLE  IF  ITS  DIMENSION  IS  DEFINED  GREATER 
THAN  5.  HOWEVER  SUBROUTINE  DRKGS  DOES  NOT  REQUIRE 
AND  CHANGE  THEM.  NEVERTHELESS  THEY  MAY  BE  USEFUL 
FOR  HANDING  RESULT  VALUES  TO  THE  MAIN  PROGRAM 
(CALLING  DRKGS)  WHICH  ARE  OBTAINED  BY  SPECIAL 
MANIPULATIONS  WITH  OUTPUT  DATA  IN  SUBROUTINE  OUTP. 

Y  -  DOUBLE  PRECISION  INPUT  VECTOR  OF  INITIAL  VALUES 
(DESTROYED). LATERON  Y  IS  THE  RESULTING  VECTOR  OF 
DEPENDENT  VARIABLES  COMPUTED  AT  INTERMEDIATE 
POINTS  X. 

DERY  -  DOUBLE  PRECISION  INPUT  VECTOR  OF  ERROR  WEIGHTS 
(DESTROYED).  THE  SUM  OF  ITS  COMPONENTS  MUST  BE 
EQUAL  TO  1.  LATERON  DERY  IS  THE  VECTOR  OF 
DERIVATIVES,  WHICH  BELONG  TO  FUNCTION  VALUES  Y  AT 
INTERMEDIATE  POINTS  X. 
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NDIM  -  AN  INPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 
EQUATIONS  IN  THE  SYSTEM. 

IHLF  -  AN  OUTPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 

BISECTIONS  OF  THE  INITIAL  INCREMENT.  IF  IHLF  GETS 
GREATER  THAN  10,  SUBROUTINE  DRKGS  RETURNS  WITH 
ERROR  MESSAGE  IHLF=1 1  INTO  MAIN  PROGRAM.  ERROR 
MESSAGE  I HL F= 1 2  OR  IHLF=13  APPEARS  IN  CASE 
P RMT ( 3 ) =  0  OR  IN  CASE  SIGN(PRMT(3) ) .NE .SIGN (PRMT(2) 
PRMT(l))  RESPECTIVELY. 

FCT  -  THE  NAME  OF  AN  EXTERNAL  SUBROUTINE  USED.  THIS 

SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDES  DERY  OF 
THE  SYSTEM  TO  GIVEN  VALUES  X  AND  Y.  ITS  PARAMETER 
LIST  MUST  BE  X,Y ,DERY.  SUBROUTINE  FCT  SHOULD 
NOT  DESTROY  X  AND  Y. 

OUTP  -  THE  NAME  OF  AN  EXTERNAL  OUTPUT  SUBROUTINE  USED. 

ITS  PARAMETER  LIST  MUST  BE  X,Y ,DERY , IHLF , NDIM, PRMT 
NONE  OF  THESE  PARAMETERS  (EXCEPT,  IF  NECESSARY, 

P RMT ( 4 ) ,PRMT(5) , . . . )  SHOULD  BE  CHANGED  BY 
SUBROUTINE  OUTP.  IF  PRMT(5)  IS  CHANGED  TO  NON-ZERO 
SUBROUTINE  DRKGS  IS  TERMINATED. 

AUX  -  DOUBLE  PRECISION  AUXIALIARY  STORAGE  ARRAY  WITH  8 
ROWS  AND  NDIM  COLUMNS. 

REMARKS 

THE  PROCEDURE  TERMINATES  AND  RETURNS  TO  CALLING  PROGRAM, IF 

(1)  MORE  THAN  10  BISECTIONS  OF  THE  INITIAL  INCREMENT  ARE 
NECESSARY  TO  GET  SATISFACTORY  ACCURACY  (ERROR  MESSAGE 
I HL F= 1 1 ) , 

(2)  INITIAL  INCREMENT  IS  EQUAL  TO  0  OR  HAS  WRONG  SIGN 
(ERROR  MESSAGES  IHLF=12  OR  IHLF=13), 

(3)  THE  WHOLE  INTEGRATION  INTERVAL  IS  WORKED  THROUGH, 

(4)  SUBROUTINE  OUTP  HAS  CHANGED  PRMT(5)  TO  NON-ZERO. 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
THE  EXTERNAL  SUBROUTINES  FCT(X ,Y ,DERY)  AND 
OUTP(X,Y, DERY, IHLF, NDIM, PRMT)  MUST  BE  FURNISHED  BY  THE  USER 

METHOD 

EVALUATION  IS  DONE  BY  MEANS  OF  FOURTH  ORDER  RUNGE-KUTTA 
FORMULAE  IN  THE  MODIFICATION  DUE  TO  GILL.  ACCURACY  IS 
TESTED  COMPARING  THE  RESULTS  OF  THE  PROCEDURE  WITH  SINGLE 
AND  DOUBLE  INCREMENT. 

SUBROUTINE  DRKGS  AUTOMATICALLY  ADJUSTS  THE  INCREMENT  DURING 
THE  WHOLE  COMPUTATION  BY  HALVING  OR  DOUBLING.  IF  MORE  THAN 
10  BISECTIONS  OF  THE  INCREMENT  ARE  NECESSARY  TO  GET 
SATISFACTORY  ACCURACY,  THE  SUBROUTINE  RETURNS  WITH 
ERROR  MESSAGE  IHLF  =  11  INTO  MAIN  PROGRAM. 

TO  GET  FULL  FLEXIBILITY  IN  OUTPUT,  AN  OUTPUT  SUBROUTINE 
MUST  BE  FURNISHED  BY  THE  USER. 

FOR  REFERENCE,  SEE 
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RALSTON/WILF,  MATHEMATICAL  METHODS  FOR  DIGITAL  COMPUTERS, 
WILEY,  NEW  YORK/LONDON,  1960,  PP.  110-120. 

SUBROUTINE  DRKGS(PRMT,Y,DERY,NDIM,IHLF,AUX) 

DIMENSION  Y(1 ) ,DERY( 1 ) ,AUX ( 8,1 ) ,A(4) ,B(4) ,C(4) ,PRMT(1) 

DOUBLE  PRECISION  PRMT.Y ,DERY,AUX,A,B,C,X,XEND,H ,AJ ,B0 ,CJ ,R1  ,R2 , 
1DELT 

DO  1  1  =  1  ,NDIM 

1  AUX(8 ,1 )  =  ,066666666666666667D0*DERY(  I ) 

X= P  RMT ( 1 ) 

XEND=PRMT (2) 

H=  P  RMT  ( 3 ) 

PRMT(5)=0.D0 
CALL  FCT(X,Y,DERY) 

ERROR  TEST 

I F( H* (XEND-X) )38,37 ,2 

PREPARATIONS  FOR  RUNGE-KUTTA  METHOD 

2  A(1 )=.5D0 

A(2)=. 29289321 881345248D0 
A(3)=l . 7071 06781 1865475D0 
A( 4)  = .  1 666666666666666700 
B ( 1 ) =2 . DO 
B ( 2 ) =1 .DO 
B  ( 3) =1 .DO 
B(4)=2.D0 
C(1  )  =  .5D0 

C (2) =. 29289321 881 345248D0 
C(3)=l .7071067811 865475D0 
C(4)=.5D0 

PREPARATIONS  FOR  FIRST  RUNGE-KUTTA  STEP 
DO  3  1  =  1  ,NDIM 
AUX ( 1  ,I)=Y(I) 

AUX ( 2,1 ) =DERY ( I ) 

AUX(3,I)=0.D0 

3  AUX(6,I)=0.D0 
IREC=0 
H=H+H 
IHLF=-1 
ISTEP=0 
IEND=0 

START  OF  A  RUNGE-KUTTA  STEP 

4  IF( (X+H-XEND)*H)7,6,5 

5  H=XEND-X 

6  IEND=1 

RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 

7  CALL  OUTP(X,Y,DERY,IREC,NDIM,PRMT) 


(XUA.UHI.f  lGa,Yaac  Y,TM>?q)23)tfla  3HITUQH8U2 
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I F  ( P  RMT ( 5 ) ) 40 ,8,40 

8  ITEST=0 

9  ISTEP=ISTEP+1 

START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
J=1 

10  AJ=A( J ) 

BJ=B(J) 

CO=C(J) 

DO  11  1=1  ,NDIM 
R1=H*DERY(I) 

R2=AJ*(R1-BJ*AUX(6 ,1 )) 

Y(I)=Y(I)+R2 

R2=R2+R2+R2 

11  AUX(6,I)=AUX(6,I)+R2-CJ*R1 
I  F(  J-4)  1 2  ,15,15 

12  J=J+1 
IF(J-3)13, 14,13 

13  X=X+.5D0*H 

14  CALL  FCT(X,Y,DERY) 

GO  TO  10 

END  OF  INNERMOST  RUNGE-KUTTA  LOOP 

TEST  OF  ACCURACY 

15  I F ( ITEST ) 1 6 , 1 6 ,20 

IN  CASE  ITEST=0  THERE  IS  NO  POSSIBILITY  FOR  TESTING  OF  ACCURACY 

16  DO  17  1  =  1  ,NDIM 

17  AUX ( 4,1) =Y  ( I ) 

ITEST=1 

I STEP= I STEP+I STEP-2 

18  I HLF=I HLF+1 
X=X-H 
H=.5D0*H 

DO  19  1  =  1  ,NDIM 
Y ( I )=AUX(1 ,1) 

DERY ( I ) =AUX(2 ,1 ) 

19  AUX ( 6 , l)=AUX(3 ,1 ) 

GO  TO  9 

IN  CASE  ITEST  = 1  TESTING  OF  ACCURACY  IS  POSSIBLE 

20  I MOD= I  STEP/ 2 

I F ( I ST  E  P - 1 MO  D - 1 MO  D ) 2 1 ,23,21 

21  CALL  FCT(X,Y,DERY) 

DO  22  1  =  1  ,NDIM 
AUX ( 5,I)=Y(I) 

22  AUX(7 ,1 )=DERY ( I) 

GO  TO  9 


. 
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COMPUTATION  OF  TEST  VALUE  DELT 

23  DELT=O.DO 

DO  24  1=1 ,  NDIM 

24  DELT=DELT+AUX(8,I)*DABS(AUX(4,I)-Y(I)) 
IF(DELT-PRMT(4) ) 28, 28, 2 5 

ERROR  IS  TOO  GREAT 

25  I F( IHLF-1 0)26 ,36 ,36 

26  DO  27  1  =  1  ,NDIM 

27  AUX(4 ,1 )=AUX(5 ,1 ) 

I STEP= I STEP+ I STEP-4 

X=X-H 

IEND=0 

GO  TO  18 

RESULT  VALUES  ARE  GOOD 

28  CALL  FCT(X,Y,DERY) 

DO  29  1=1  ,NDIM 
AUX ( 1  ,1 )=Y( I ) 

AUX(2 ,1)=DERY(  I) 

AUX(3,I)=AUX(6,I) 

Y  ( I )=AUX(5 ,1 ) 

29  DERY ( I )=AUX(7 ,1 ) 

CALL  OUTP ( X-H  ,Y ,DERY , IHLF, NDIM, PRMT) 

I F ( P  RMT ( 5 ) )40 ,30 ,40 

30  DO  31  1  =  1  ,NDIM 
Y ( I ) =AUX ( 1 ,1) 

31  DERY ( I )=AUX(2 ,1 ) 

I REC= I HLF 

I F( I  END ) 32 ,32 ,39 

INCREMENT  GETS  DOUBLED 

32  IHLF=I HLF-1 

I STEP= I STEP/2 
H=H+H 

I F( IHLF) 4 ,33 ,33 

33  I MOD= I  STEP/2 

I F( ISTEP-IM0D-IM0D)4 ,34 ,4 

34  IF(DELT- ,02D0*PRMT(4) ) 35 ,35 ,4 

35  I HLF=I HLF-1 
ISTEP=ISTEP/2 
H=H+H 

GO  TO  4 

RETURNS  TO  CALLING  PROGRAM 

36  I HL F= 1 1 

CALL  FCT ( X , Y ,DERY) 

GO  TO  39 

37  I HL F= 1 2 
GO  TO  39 

38  I HL F=1 3 

39  CALL  OUTP(X,Y, DERY, IHLF, NDIM, PRMT) 

40  RETURN 
END 
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